МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ 
им. М.В. ЛОМОНОСОВА
-----------------------------------------------------------------------------------------------------------------------
Научно-исследовательский вычислительный центр

О.Б. Арушанян, Н.И. Волченскова

РЕШЕНИЕ ОБОБЩЕННОЙ АЛГЕБРАИЧЕСКОЙ
ПРОБЛЕМЫ СОБСТВЕННЫХ ЗНАЧЕНИЙ
ДЛЯ СИММЕТРИЧНЫХ И ЭРМИТОВЫХ МАТРИЦ
НА РАСПРЕДЕЛЕННОЙ ПАМЯТИ

Учебное пособие

Москва, 2017

О.Б. Арушанян, Н.И. Волченскова

Решение обобщенной алгебраической проблемы собственных значений для симметричных и эрмитовых матриц на распределенной памяти

(Учебное пособие)

Оглавление

1. Введение                                                                                                                               4
    2. Основные правила решения обобщенной проблемы собственных значений 
       для симметричных и эрмитовых матриц с помощью комплекса PARALG                   5
2.1.	Общее описание организации и структуры комплекса                                        5
2.2.	Применение пакета BLACS для организации параллельных процессов           7
2.3.	Задание параметров при обращении к целевым программам комплекса           8
2.3.1. Глобальные и локальные объекты и параметры параллельных программ        8
2.3.2. Дескрипторы глобальных массивов                                                                     10
2.4.	Алгоритмы распределения матриц по решетке процессов                                 11
2.4.1. Блочно-циклическое отображение плотных матриц в локальную 
          память параллельных процессов                                                                           11
2.4.2. Пример блочно-циклического распределения плотной матрицы по 
          решетке процессов                                                                                                  14
2.5.	Алгоритмы решения обобщенной проблемы собственных значений 
              для симметричных и эрмитовых матриц                                                               16
3.	Практические сведения по использованию параллельных программ комплекса         18
3.1.	Предварительные действия, необходимые при обращении к 
                программам комплекса                                                                                 18
3.2.	Программирование распределения матриц по решетке процессов         21
3.3.	О документировании и примерах использования 
                 параллельных программ                                                                              22
3.4.	Пример использования программы решения обобщенной 
                 проблемы собственных значений для симметричных матриц                26
3.5.	Общий список программ комплекса PARALG для решения 
                 обобщенной проблемы собственных значений с симметричными
                 и эрмитовыми матрицами                                                                           28
3.6.	Запуск программ комплекса в ОС Linux и получение результатов
                на суперкомпьютере «Чебышев»                                                                29
3.7.	Запуск программ комплекса в ОС Linux и получение результатов
                 на суперкомпьютере «Ломоносов-1»                                                         31
Приложение 1. Описание подпрограммы для вычисления всех собственных 
                            значений и собственных векторов вещественной 
                            симметричной матрицы  в обобщенной проблеме                                   33 
Приложение 2. Головная программа для вызова подпрограммы из Приложения 1
                          при решении обобщенной проблемы собственных значений                   42

Литература                                                                                                                     47

1. Введение
Уровень и сложность решаемых в настоящее время научных и производственных задач таков, что, как правило, требуется использование высокопроизводительных суперкомпьютеров, в частности, с распределенной памятью. Подготовка специалистов в этой области в МГУ им. М.В.Ломоносова ведется на базе Суперкомпьютерного Центра (СКЦ МГУ), состоящего в настоящее время из трех высокопроизводительных суперкомпьютеров: СКИФ МГУ «Чебышев»(60 Тфлоп/с), «Ломоносов-1» (1,7 Пфлоп/с) и «Ломоносов-2» (2,96 Пфлоп/с).
Описываемый в настоящем учебном пособии комплекс программ в основном был реализован с использованием суперкомпьютера СКИФ МГУ «Чебышев», а затем перенесен на суперкомпьютер Ломоносов-1.
Проблема собственных значений является одной из наиболее часто встречающихся численных задач. В случае матриц больших размеров требуются параллельные алгоритмы, позволяющие использовать многопроцессорные вычислительные системы.  
В настоящем пособии рассматривается частный случай (один из видов указанной задачи) - обобщенная проблема собственных значений для симметричных и эрмитовых матриц.
Математическая запись такой задачи состоит в следующем.  
Пусть A и B являются вещественными симметричными или эрмитовыми матрицами порядка n, а матрица B является также положительно определенной.
    Обобщенная симметричная проблема собственных значений состоит в нахождении собственных значений λ и соответствующих собственных векторов z, удовлетворяющих одному из следующих видов уравнений:

     (1) 				A z  = λ B z ,

     (2)				A B z  = λ z ,

     (3)				B A z  = λ z .

При этом для всех видов уравнений (1), (2) или (3) искомые значения λ всегда являются вещественными.

 Высокоэффективное (в смысле производительности) решение таких задач может быть выполнено с использованием специализированных пакетов программ, разработанных высококвалифицированными специалистами в этой области. Одним из таких наиболее широко известных пакетов является пакет ScaLAPACK (Scalable Linear Algebra PACKage) [15,16], позволяющий решать задачи с матрицами больших размеров на распределенной памяти с использованием технологии MPI(Massage Passing Interface) [1,2].
В НИВЦ МГУ создан русскоязычный аналог этого пакета, получивший название «Комплекс параллельных программ PARALG». В его основе лежат те же принципы реализации параллельных алгоритмов, которые приняты в ScaLAPACKе, и при реализации целевых программ комплекса PARALG использовались базовые модули пакета ScaLAPACK. В комплексе PARALG содержатся и задокументированы (т.е.  снабжены инструкциями по использованию) именно готовые целевые программы решения задач линейной алгебры (в том числе для решения обобщенной проблемы собственных значений). Комплекс содержит большее количество целевых программ, чем в пакете ScaLAPACK; полный список их приводится в п.3.5 настоящего учебного пособия. Базовые же модули (как составные части целевых программ) вызываются непосредственно из библиотеки ScaLAPACKа.
Общие правила и принципы, приводимые в теоретической части предлагаемого учебного пособия (п.2), в одинаковой мере справедливы как для пакета ScaLAPACK, так и для комплекса PARALG. Использование готовых целевых программ комплекса PARALG значительно облегчает эффективное решение указанных задач.  
Тем не менее, у неспециалистов в данной области (как показывает опыт) это может вызвать определенные трудности, так как требует знакомства с недостаточно известными действиями и понятиями пакета ScaLAPACK, связанными с обеспечением эффективности параллельных вычислений. Сюда относятся такие, например, особенности параллельной реализации алгоритмов как:
- организация параллельных процессов в так называемую (виртуальную) решетку процессов и сопоставление с выбранной решеткой определенного численного значения параметра, который однозначно характеризует выбранную решетку и называется указателем контекста или просто контекстом;
- специальное распределение блоков матриц по параллельным процессам решетки;
- понятие и правильное определение так называемых дескрипторов матриц и связанное с этим представление о локальных и глобальных параметрах целевых программ комплекса PARALG.
Кроме того, требуется овладеть правильным использованием подпрограмм пакета BLACS (Basic Linear Algebra Communication Subprograms) [19, 20], входящего в состав пакета ScaLAPACK и выполняющего на более высоком уровне функции примитивов MPI. Указанные особенности обсуждаются в п.2 данного учебного пособия. Кроме того, в п.2 приводится краткое описание математических аспектов решения обобщенной проблемы собственных значений для симметричных и эрмитовых матриц.
В п.3. разбирается пример программирования при решении обобщенной проблемы собственных значений для симметричных и эрмитовых матриц с помощью программ комплекса PARALG. Кроме того, описываются
- запуск соответствующих задач в ОС Linux на суперкомпьютере СКИФ МГУ «Чебышев»;
- запуск соответствующих задач в ОС Linux на суперкомпьютере «Ломоносов-1».
2.	Основные правила решения обобщенной проблемы собственных значений для симметричных и эрмитовых матриц с помощью комплекса PARALG 

2.1.Общее описание организации и структуры комплекса

            Комплекс программ PARALG относится к разряду “параллельных предметных библиотек”, при обращении к которым пользователю не приходится явно применять какие-либо конструкции специальных языков параллельного программирования или интерфейсов, поддерживающих взаимодействие параллельных процессов. Все необходимые действия по совместной работе параллельных процессов при решении задач линейной алгебры выполняются подпрограммами из пакетов PBLAS (Parallel Basic Linear Algebra Subprogramms)[17, 18] и BLACS (Basic Linear Algebra Communication Subprogramms), являющихся составными частями пакета ScaLAPACK.

     Пакет BLACS является набором подпрограмм, предназначенных для использования методов распараллеливания при решении задач линейной алгебры. Этот пакет разработан для компьютеров с распределенной памятью и обеспечивает запуск параллельных процессов и их взаимодействие посредством обмена сообщениями. При этом пакет BLACS освобождает пользователей от изучения предназначенных для этих целей комплексов стандартных примитивов, включенных в MPI.  

Работа подпрограмм пакета BLACS опирается на представление о так называемой «решетке процессов» (process grid).

Предполагается, что параллельные процессы организованы в двумерную (или одномерную) решетку процессов, в которой каждый из процессов идентифицируется своими координатами в решетке и хранит определенные части обрабатываемых матриц и векторов. Для работы с параллельными процессами в состав пакета BLACS включены подпрограммы, реализующие следующие основные функции:
- построение, изменение и получение сведений о решетке процессов;
- обмен сообщениями (передача матриц или их частей) между двумя процессами;
- передача сообщений от одного процесса сразу многим процессам;
- выполнение необходимых итоговых действий после получения параллельными процессами частичных результатов (например, вычисление итоговых сумм, максимумов или минимумов по всем процессам).
Другим пакетом, к которому обращаются программы комплекса PARALG, является пакет PBLAS, содержащий версии для распределенной памяти подпрограмм (первого, второго и третьего уровней) пакета BLAS (Basic Linear Algebra Subprogram) [22, 23]. В этом пакете содержатся версии подпрограмм, реализующих базовые операции линейной алгебры (например, действия над векторами и скалярами, умножение матрицы на вектор или перемножение двух матриц и др.).
Разработчики пакета PBLAS реализовали его таким образом, чтобы интерфейс его подпрограмм был максимально похожим на интерфейс подпрограмм пакета BLAS.  
Таким образом, целевые программы описываемого комплекса PARALG (т.е. программы, которые имеют самостоятельное значение при решении прикладных задач) могут содержать обращения к подпрограммам указанных выше пакетов BLACS, PBLAS и BLAS, а также к базовым и вспомогательным подпрограммам пакета ScaLAPACK.
К базовым подпрограммам пакета ScaLAPACK относятся подпрограммы, которые не используются независимо от других программ комплекса PARALG. Они всегда выступают в роли вызываемых из других подпрограмм. Например, подпрограмма умножения матрицы общего вида на ортогональную матрицу, полученную в результате QR - разложения другой подпрограммой.
К вспомогательным подпрограммам относятся подпрограммы, называемые tools routines. Эти подпрограммы выполняют различные вспомогательные функции, например: 
- выдача параметров машинной арифметики конкретной ЭВМ;
- выдача диагностических сообщений в стандартной форме;
- подсчет количества строк или столбцов распределенной по параллельным процессам матрицы, расположенных на одном из этих процессов и др.
Целевые программы комплекса PARALG решают следующие задачи из разделов линейной алгебры:
- решение систем линейных алгебраических уравнений, 
- решение проблемы собственных значений (линейной и обобщённой),
- сингулярное разложение матриц. 
Настоящее пособие предназначено для обучения практическим навыкам использования программ комплекса для решения обобщенной проблемы собственных значений для симметричных и эрмитовых матриц.
2.2.	Применение пакета BLACS для организации параллельных процессов

Подпрограммы пакета BLACS разработаны специально для обеспечения эффективного решения на параллельных компьютерах с распределенной памятью задач линейной алгебры, в которых основным обрабатываемым объектом является матрица.  
При распределении таких двумерных массивов по параллельным процессам удобнее мысленно организовать параллельные процессы в виде логической двумерной решетки процессов. Тогда, например, все части строки распределенной матрицы могут быть найдены в процессах, относящихся к строке решетки процессов, а части столбца - в процессах, относящихся к столбцу решетки процессов. 
В принятой в BLACS логической структуре параллельных процессов каждый процесс идентифицируется с помощью пары чисел, являющихся координатами процесса в решетке процессов.  
На рисунке ниже изображена решетка из 8 процессов, состоящая из двух строк и четырех столбцов (2 x 4):
                    0                       1                     2                     3
------------------------------------------------------------------------------
0   |                         |                         |                        |                      |
------------------------------------------------------------------------------
1   |                         |                         |                        |                      |
------------------------------------------------------------------------------

Этим процессам ставятся в соответствие координаты решетки:

( 0,0 ), ( 0,1 ), ( 0,2 ), ( 0,3 ), ( 1,0 ), ( 1,1 ), ( 1,2 ), ( 1,3 )

      В общем случае, если решетка процессов имеет R строк и C столбцов, то любой процесс в решетке обозначается парой (i, j ), где 0 ≤ i < R и 0≤ j < C. Общее число процессов равно R*C.

      В компьютерах с распределенной памятью каждый параллельный процесс работает со своей локальной памятью, а обмен информацией между процессами осуществляется посредством передачи сообщений через связывающую их сеть.

     Пакет BLACS, предназначенный для таких компьютеров, опирается на использование
широко известных и распространненых интерфейсов обмена сообщениями MPI (Massage
Passing Interface) или PVM (Parallel Virtual Machine). Они содержат не только операции обмена между одним посылающим и одним принимающим сообщение процессом, но также и массовые операции обмена, когда, например, один процесс может рассылать сообщение сразу всем другим параллельным процессам, и, наоборот, один из процессов может получать сообщения от всех других.

В пакете BLACS такие операции называются контекстными операциями (scope operations), а про процессы, принимающие в них участие, говорят, что они находятся внутри операционного контекста (operation's scope). Организация процессов в виде логической решетки, позволяет пакету BLACS естественным образом расширить виды контекстов, а именно, объявлять (задавать) в качестве контекста, т.е. области действия массовых операций обмена сообщениями, не только все параллельные процессы, но также и все процессы, принадлежащие только некоторой строке решетки процессов или некоторому ее столбцу. Это обеспечивает программисту в области линейной алгебры более удобные и гибкие возможности распараллеливания.

     Использование понятия решеток процессов удобно не только для программирования, но и позволяет повысить эффективность (скорость) вычислительного процесса. Однако это произойдет только в том случае, если в базовой машине параллельные процессы связаны между собой с помощью, по крайней мере, двумерных сетей.

В большинстве современных суперкомпьютеров это условие выполнено, но при использовании сети Ethernet такое преимущество недостижимо.

2.3.	Задание параметров при обращении к целевым программам комплекса

2.3.1.	Глобальные и локальные объекты и параметры программ

     Программы комплекса PARALG составлены на языке ФОРТРАН - 77 в стиле SMPD (Single Program Multiple Data). Это означает, что одна и та же исполняемая программа загружается во все параллельные процессы, заказанные для решения задачи в команде запуска на счет. В то же время, каждый процесс занимается обработкой своих данных, составляющих некоторую часть общих данных задачи. Эти данные отдельного процесса (их называют локальными данными) непосредственно недоступны другим процессам и хранятся в отдельной (локальной) памяти каждого процесса.  Необходимый обмен данными между разными процессами производится только посредством передачи сообщений (технология MPI [1, 2, 8 - 11]).
  Из сказанного следует, что исходную матрицу  A не требуется хранить целиком в оперативной памяти процессора как единый массив. На практике такие матрицы (векторы) в виде единого целого могут храниться во внешних файлах (на дисковой памяти). По этой причине такую матрицу  A называют глобальной матрицей.
  Поскольку перед началом работы целевых параллельных программ комплекса исходные глобальные матрицы (или векторы) разделяются на части и распределяются по параллельным процессам, то в дальнейшем нередко употребляется термин «распределенная глобальная матрица (вектор)». Часть глобальной матрицы (вектора), представленная в локальной памяти одного из параллельных процессов, называется локальной частью матрицы (вектора). Совокупность всех локальных частей на всех, используемых при решении задачи параллельных процессах, и составляет распределенную глобальную матрицу (вектор).
  Поскольку каждый процесс работает со своей локальной частью данных (в нашем случае фрагментом глобальной матрицы), то в головной фортранной программе достаточно выделить (заказать) память, необходимую только для максимальной локальной части этой матрицы из всех локальных частей, обрабатываемых на параллельных процессах.
Локальные части, обрабатываемые на каждом из процессов, могут быть как одинакового, так и разного размера (см. подробнее пп. 2.4 и 3.2). Для того чтобы каждый из процессов мог разместить свою локальную часть данных (матрицы), заказываемая область локальной памяти должна быть не меньше максимальной локальной части данных.
Таким образом, в головной фортранной программе, обращающейся к целевой программе комплекса, не требуется резервировать память для всего глобального объекта (матрицы или вектора) целиком, а только для их локальных частей.
Целиком эти объекты существуют, как правило, во внешних файлах. Для сохранения информации о конкретном способе разбиения глобального объекта на блоки и размещении блоков в локальной памяти всех процессов вводится объект, называемый дескриптором глобального массива. Дескриптор представляет собой одномерный массив, состоящий из 9 или 7 элементов (в зависимости от типа дескриптора). В дескрипторе хранится информация о размерах глобального массива (число строк и столбцов матрицы), о размерах блоков, на которые она разбивается, о номере процесса, в память которого распределяется первый блок глобальной матрицы (левый верхний угол), а также некоторая другая информация. Подробнее о дескрипторах см.  в п. 2.3.2.
Учитывая общие для комплекса правила распределения блоков матриц (векторов) по параллельным процессам, по информации в дескрипторах всегда можно определить, в локальной памяти какого процесса и где именно находится элемент глобальной матрицы с глобальными индексами  i  и  j (см. подробнее в п. 2.4.1).  
Непосредственное размещение исходных данных в локальной памяти каждого процесса должен выполнять сам пользователь (в головной программе перед обращением к целевой программе комплекса) в соответствии с установленными правилами для матриц разных видов: плотных общего вида или симметричных, ленточных или трехдиагональных.
Организация решетки процессов, обмен информацией между процессами, а также другие действия, связанные с процессами, выполняются с помощью подпрограмм пакета BLACS. Пользователю, однако, в большинстве случаев не требуется знание всех возможностей этого пакета. При решении своей задачи он должен сделать в головной программе лишь несколько стандартных обращений к нескольким подпрограммам пакета BLACS.
Наглядные примеры этого можно найти, например, в предлагаемом пользователям описании целевой параллельной подпрограммы комплекса (в разделе «Пример использования»; см. Приложение 1), а также в фортранном тексте тестового примера к этой подпрограмме (см. Приложение 2). Более подробно описание решеток процессов и правил работы с ними пользователей см. в пп. 2.3.2, 2.4 и 3.1.  
Из сказанного выше следует, что формальные параметры целевых программ комплекса могут быть глобальными или локальными.
Глобальными называются параметры подпрограмм, которые относятся к указанным выше глобальным объектам. Например, глобальными являются число строк M или столбцов N исходной матрицы  A или число наддиагоналей BWU в матрице A, если она является ленточной.
Локальными называются параметры, которые характеризуют объекты, определяемые в локальной памяти отдельного процесса. Например, указатель на локальную память, занимаемую локальной частью глобальной матрицы  A или глобального вектора  B правых частей линейной системы уравнений.
Упоминавшиеся выше дескрипторы массивов характеризуются одновременно и как глобальные, и как локальные параметры подпрограмм, поскольку среди их элементов, определяющих в основном характеристики глобальных массивов, есть характеристика, описывающая размещение части глобального массива в локальной памяти процесса (ведущая размерность локального массива). 

2.3.2.	Дескрипторы глобальных массивов

Как уже упоминалось выше, для каждого глобального массива (матрицы или вектора), который необходимо разбить на блоки (с последующим размещением каждого из них в локальной памяти параллельных процессов, участвующих в решении задачи), вводится специальный объект, называемый дескриптором массива. Дескриптор представляет собой одномерный массив, состоящий из 9 или 7 элементов целого типа (в смысле языка ФОРТРАН). Он предназначен для хранения информации, конкретизирующей способ разбиения глобального массива на блоки и их распределение по всем параллельным процессам.
Дескрипторы являются обязательными параметрами целевых программ комплекса.  Дескрипторы принято обозначать идентификатором DESC с добавлением в конце имени массива (матрицы или вектора). Например, DESCA означает дескриптор матрицы  A и DESCB - дескриптор вектора  B.
Дескрипторы бывают трех типов, которые обозначаются целыми числами: 1, 501, 502.
Для плотных глобальных матриц  A, которые распределяются по двумерной решетке процессов, используется дескриптор, имеющий тип 1 и состоящий из 9 элементов, за каждым из которых закреплено свое символическое имя (см. таблицу ниже).  В таблице указываются номер (индекс) каждого элемента в массиве, его символическое имя (которое оканчивается символом подчеркивания «_»), а также его смысл. Идентификатор, стоящий в символическом имени элемента после символа подчеркивания «_», является идентификатором глобального массива, описываемого данным дескриптором. Например, N_A обозначает элемент дескриптора, который содержит число столбцов глобальной матрицы  A.

Таблица элементов дескриптора плотной глобальной матрицы A.

N           Имя                                                   Смысл
         1.      DTYPE_A  -  тип дескриптора (для плотной матрицы  A  DTYPE_A = 1)

2.      CTXT_A    - обозначение контекста BLACS'а, соответствующего выбранной
                               решетке процессов, по которой распределяется глобальная матрица  
                              A; целое значение, устанавливаемое подпрограммой из пакета 
                              BLACS; см. п. 3.1.

3.       M_A         - число строк в глобальной матрице  A

4.       N_A         - число столбцов в глобальной матрице  A

5.       MB_A      - число строк в блоках, на которые разбивается глобальная матрица A

6.       NB_A       - число столбцов в блоках, на которые разбивается глобальная матрица 
                           A

7.       RSRC_A  - номер строки процесса в решетке процессов, куда был распределен
                             первый элемент глобальной матрицы  A (левый верхний угол); как 
                             правило удобнее всего распределять первый элемент в процесс с 
                             номером (0,0)

8.      CSRC_A   - номер столбца процесса в решетке процессов, куда был распределен 
                             первый элемент глобальной матрицы  A (левый верхний угол); как 
                             правило, удобнее всего распределять первый элемент в процесс с 
                             номером (0,0)

9.     LLD_A      - ведущая размерность локального массива: 
                             LLD_A ≥ MAX(1, LOCr(M_A))

В этой таблице и далее число строк локального массива обозначается LOCr(M_A). 
Заметим, что записи DESCA(3) и DESCA(M_A) означают один и тот же элемент дескриптора матрицы A, содержащий число строк этой матрицы.

2.4.	Алгоритмы распределения матриц по решетке процессов

2.4.1.	Блочно-циклическое отображение плотных матриц в локальную память параллельных процессов

Как уже говорилось выше, перед обращением к программам комплекса PARALG пользователю необходимо разделить исходную матрицу на блоки и распределить ее по параллельным процессам. Для матриц разных видов приняты свои способы разбиения и распределения матриц. Здесь мы рассмотрим наиболее общий случай прямоугольной плотной матрицы общего вида 
Для равномерной загрузки параллельных процессов и эффективного использования подпрограмм пакета BLAS при выполнении операций над плотными глобальными матрицами были разработаны специальные способы разбиения таких матриц на прямоугольные блоки MB x NB и распределения этих блоков по двумерной решетке параллельных процессов. Простой пример конкретного разбиения плотной матрицы на блоки и распределения их по двумерной решетке процессов, называемого блочно-циклическим отображением, приводится в п. 2.4.2.
Здесь же мы рассмотрим, каким образом блоки, отображенные на один и тот же процесс, располагаются и хранятся в локальной памяти этого процесса. Другими словами, мы выпишем точные формулы, которые связывают элемент глобальной матрицы, определяемый глобальными индексами I и J (т.е. A (I, J)), с координатами процесса, владеющего этим элементом, в решетке процессов (Pr, Pc) и с расположением этого элемента матрицы в локальной памяти этого процесса.  
Рассмотрим сначала, для простоты, эти связи на примере одномерного случая, т.е.  размещение одномерного глобального массива длины N, разбитого на блоки длины NB, по одномерной решетке из P процессов, перенумерованных начиная с 0 до (P - 1).  Элементы самого глобального массива нумеруются от 1 до N. Прежде всего, массив делится на смежные блоки размера NB. Когда N не делится на NB нацело, последний блок массива элементов будет содержать только mod (N, NB) элементов вместо NB элементов. Блоки, на которые разбивается исходный массив, нумеруются (как и процессы) начиная с 0. Они распределяются по процессам подобно тому, как сдается колода карт участникам игры - по кругу (т.е. циклически).  
Другими словами, мы предполагаем, что процесс с номером 0 получает первый блок (с номером 0), k-й блок связывается с процессом, имеющим координату (номер) mod (k, P). Блоки, связанные с одним и тем же процессом, хранятся в памяти рядом (в смежных, прилегающих частях). Отображение элемента массива с глобальным индексом I определяется следующим соотношением:
I   =   k * NB + x   =   ( m * P + p ) * NB + x  ,
     где

I - глобальный индекс элемента глобального массива; 
m - локальная координата блока, в котором этот элемент расположен; 
р - координата процесса, владеющего этим блоком; 
х - локальная координата элемента внутри того блока, где находится элемент 
     глобального массива с индексом  I.
Легко установить следующие соотношения между этими переменными:
p   =   [ ( I - 1) / NB ] mod P ,
m  =   [ ( I - 1) / ( P * NB ) ] ,
x   =   mod ( I - 1, NB ) + 1 .

Эти уравнения позволяют определить локальную информацию (т.е. локальный индекс m * NB + x), так же как и координату процесса  р, соответствующую глобальному элементу, идентифицируемому его глобальным индексом  I, и наоборот.  В  таблице ниже показано отображение в локальную память при разбиении массива на блоки, когда  P = 2 ,  N = 16 и  NB = 8.
          I                   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
p                  0   0   0   0   0   0   0   0   1   1    1    1    1    1    1    1
m                 0   0   0   0   0   0   0   0   0   0    0    0    0    0    0    0
x                  1   2   3   4   5   6   7   8   1   2    3    4    5    6    7    8
m*NB + x   1   2   3   4   5   6   7   8   1   2    3    4    5    6    7    8

Вовсе не обязательно всегда делать распределение блоков начиная с процесса с номером 0. Иногда бывает полезно начать распределение данных с процесса, имеющего отличную от нуля координату  SRC. В этом случае соотношения становятся такими:
p   =   ( SRC + [ ( I - 1) / NB ] ) mod P ,
m  =   [ ( I - 1) / ( P * NB ) ] ,
x   =   mod ( I - 1, NB ) + 1 .
Ниже в таблице показано размещение блоков при  P = 2,  SRC = 1,  NB = 3 и  N = 16.
           I                  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
p                  1   1   0   0   0   1   1   1   0   0    0    1    1    1    0    0
m                 0   0   0   0   0   0   1   1   1   1    1    1    2    2    2    2
x                  1   2   3   1   2   3   1   2   3   1    2    3    1    2    3    1
m*NB + x   1   2   3   1   2   3   4   5   6   4    5    6    7    8    9    7

В двумерном случае предположим, что матрица разделена на блоки размера MB x NB и что первый блок передан процессу с координатами (RSRC, CSRC). Формула, приведенная выше для одномерного случая, должна использоваться повторно независимо для каждого измерения решетки процессов  Pr x Pc. Например, элемент матрицы ( I, J ) находится в процессе с координатами ( Pr, Pc) внутри локального блока ( m, n ) в позиции ( x, y ), задаваемой формулами
( m, n )    =   ( [ ( I - 1) / ( Pr * MB ) ], [ ( J - 1) / ( Pc * NB ) ] ) ,

( Pr, Pc)   =   ( ( RSRC + [ ( I - 1) / MB ] ) mod Pr, ( CSRC + [ ( J - 1) / NB ] ) mod Pc )

( x, y )    =   ( mod ( I - 1, MB ) + 1, mod ( J - 1, NB ) + 1 )

    Эти формулы показывают, как матрица  А размера M_A x N_A отображается и хранится на решетке процессов. Матрица сначала разбивается на блоки размера 
MB_A x NB_A начиная с ее верхнего левого угла. Эти блоки затем равномерно распределяются по решетке процессов циклическим образом.

    Каждый процесс владеет набором блоков, которые хранятся рядом (смежно) по столбцам в двумерном массиве, расположенном в памяти «по столбцам».

    Это соглашение о локальном размещении позволяет эффективно использовать иерархию локальной памяти посредством вызова пакета BLAS для подмассивов, которые могут быть больше, чем один блок размера MB_A x NB_A.

    На рисунках ниже представлены матрица, разделенная на блоки размером  2 x 2, и отображение матрицы  5 x 5 на решетку процессов размером 2 x 2 (т.е. M_A = N_A = 5, Pr= Pc = 2 и MB_A = NB_A = 2). Локальные элементы каждого столбца матрицы хранятся рядом в памяти каждого процесса.

a11  a12          a13  a14          a15
a21  a22          a23  a24          a25
----------------------------------------
a31  a32          a33  a34          a35
a41  a42          a43  a44          a45
----------------------------------------
a51  a52          a53  a54          a55

    Цифры 0 и 1 в левой и верхней частях приведенной ниже таблицы означают номера строк и столбцов решетки процессов соответственно. 

          0                           1
     ----------------------------------------------------
                  a11  a12    a15          a13  a14
           0     a21  a22    a25          a23  a24
------------------------------------------
a51  a52    a55          a53  a54
     ----------------------------------------------------
           1     a31  a32    a35          a33  a34
a41  a42    a45          a43  a44

    Важной задачей для пользователя является определение числа строк и столбцов плотной глобальной матрицы, которое получает каждый конкретный процесс. Для выполнения этой функции существует служебная подпрограмма NUMROC.

    Для локального числа строк используется обозначение LOCr ( ), а для локального числа столбцов обозначение LOCc ( ). Значения LOCr ( ) и LOCc ( ), получаемые подпрограммой NUMROC, являются результатом точных вычислений.

    Однако если требуется понять общую идею вычисления размера локального массива, то можно проделать следующее грубое вычисление верхней границы этой величины: 
    а) верхняя граница для  LOCr ( ) оценивается по формуле 
LOCr( )   =   [ [ M_A / MB_A] / Pr ] * MB_A

         б) величина LOCc( ) оценивается по формуле
LOCc( )   =   [ [ N_A / NB_A ] / Pc ] * NB_A

    Заметим, что эти вычисления могут привести к очень большой переоценке величины действительно требующегося пространства.

2.4.2.	Пример блочно-циклического распределения плотной матрицы по
                решетке процессов

Приводимый в настоящем разделе пример является простой иллюстрацией блочно-циклического распределения плотной глобальной матрицы  A по двумерной решетке процессов.
В соответствии с принятой схемой, сначала плотная матрица  A размера M x N разделяется на блоки размера MB x NB начиная с левого верхнего угла этой матрицы. Эти блоки затем равномерно распределяются по каждому измерению решетки процессов. Математические соотношения, соответствующие такой схеме распределения, приводятся в п. 2.4.1.
Таким образом, каждый процесс владеет набором блоков, которые расположены в его локальной памяти рядом в двумерном массиве, хранящемся по столбцам (см. ниже).  На приведенном ниже рисунке показано разделение матрицы  A размером 9 x 9 на блоки размером 2 x 2.
a11  a12          a13  a14          a15  a16          a17  a18          a19
a21  a22          a23  a24          a25  a26          a27  a28          a29
a31  a32          a33  a34          a35  a36          a37  a38          a39
a41  a42          a43  a44           a45  a46          a47  a48         a49
a51  a52          a53  a54          a55  a56          a57  a58          a59
a61  a62          a63  a64          a65  a66          a67  a68          a69
a71  a72          a73  a74          a75  a76          a77  a78          a79
a81  a82          a83  a84          a85  a86          a87  a88          a89
a91  a92          a93  a94          a95  a96          a97  a98          a99

     Далее показано отображение полученных блоков на решетку процессов размером 2 * 3 (двумерное блочно-циклическое распределение данных). Чтобы понять, как распределятся изображенные выше блоки по процессам решетки, сначала напишем на каждой из клеток (блоков) рисунка координаты процессов, куда они должны быть распределены в соответствии с установленным блочно-циклическим распределением, описанным в п. 2.4.1.

         (0,0)     (0,1)     (0,2)     (0,0)     (0,1)
         (1,0)     (1,1)     (1,2)     (1,0)     (1,1)
         (0,0)     (0,1)     (0,2)     (0,0)     (0,1)
         (1,0)     (1,1)     (1,2)     (1,0)     (1,1)
         (0,0)     (0,1)     (0,2)     (0,0)     (0,1)
     Пусть первый блок распределяется в процесс (0, 0). Тогда все блоки, расположенные с ним в той же строке, будут распределяться в ту же строку решетки процессов, т.е. первая координата в этих клетках будет равна 0.
     Вторая же координата (столбца решетки) будет циклически изменяться: 0, 1, 2, 0, 1, так как столбцов в решетке только 3 (с координатами 0, 1, 2). Все блоки, расположенные с первым блоком в том же самом столбце, будут распределяться в тот же самый столбец решетки процессов. Иными словами, вторая координата в этих клетках будет равна 0. Первая же координата (строки решетки) будет циклически изменяться: 0, 1, 0, 1, 0, поскольку число строк в решетке равно только 2 (с координатами 0 и 1).
Другими словами, каждая координата независимо от другой (в строках - слева направо, а в столбцах - сверху вниз) циклически изменяется.  
Если теперь мы соединим вместе (пристыкуем) все блоки, имеющие одинаковые координаты процесса, то получим массив элементов, который и должен расположиться в локальной памяти процесса с такими координатами.
Сборка клеток (блоков) с одинаковыми координатами делается таким образом:
- из этих клеток выбирается расположенная левее и выше всех; 
- к ней пристыковываются снизу клетки, расположенные в этом же столбце, тем самым получается первый столбец блоков; 
- затем берется клетка, расположенная правее в той же строке клеток, что и самая первая; к ней пристыковываются все расположенные ниже в том же столбце, тем самым получаем следующий столбец клеток; 
- после того, как собраны все столбцы клеток с одинаковыми координатами, пристыковываем их по-очереди справа к первому столбцу с такими же координатами.
Тем самым получаем единый двумерный массив, расположенный в локальной памяти процесса с указанными в клетках координатами.
Ниже представлено расположение элементов исходной матрицы во всех процессах решетки после завершения процесса блочно - циклического распределения. Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.
                 0                    |           1                         |            2
          --------------------------------------------------------------------------------------
a11  a12         a17  a18          a13  a14        a19          a15  a16
a21  a22         a27  a28          a23  a24        a29          a25  a26
-------------------------------------------------------------------------------
          0      a51  a52          a57  a58          a53  a54         a59          a55  a56
a61  a62         a67  a68          a63  a64         a69          a65  a66
-------------------------------------------------------------------------------
a91  a92         a97  a98          a93  a94         a99          a95  a96
          --------------------------------------------------------------------------------------
a31  a32         a37  a38          a33  a34         a39          a35  a36
a41  a42         a47  a48          a43  a44         a49          a45  a46
          1      --------------------------------------------------------------------------------------
a71  a72         a77  a78          a73  a74         a79          a75  a76
a81  a82         a87  a88          a83  a84         a89          a85  a86

Ниже в таблице указаны характеристики локальных массивов для каждого из процессов решетки.
Координаты              LLD_A       LOCr ( M_A )       LOCc ( N_A )
процесса
-----------------------------------------------------------------------------------------------------
(0,0)                      5                     5                             4
(0,1)                      5                     5                             3
(0,2)                      5                     5                             2

(1,0)                      4                     4                             4
(1,1)                      4                     4                             3
(1,2)                      4                     4                             2

Число строк LOCr и число столбцов LOCc матрицы  A, которыми владеет конкретный процесс, могут отличаться у разных процессов в решетке. Подобно этому, для каждого процесса в решетке процессов существует ведущая локальная размерность LLD. Ее величина может быть различной для разных процессов в решетке процессов.  Например, как мы можем видеть на рисунке выше, локальный массив, хранящийся в строке решетки процессов с номером 0, должен иметь ведущую локальную размерность LLD не меньше 5, а хранящийся в строке с номером 1 - не меньше 4.
2.5.	Алгоритмы решения обобщенной симметричной проблемы собственных
значений

      Пусть A и B являются вещественными симметричными или эрмитовыми матрицами порядка  n, а матрица B является также и положительно определенной.  Обобщенная симметричная проблема собственных значений состоит в нахождении собственных значений λ  и соответствующих собственных векторов  z, удовлетворяющих одному из видов уравнений

      (1)           A z   =   λ B z ,

      (2)           A B z  =  λ  z ,

      (3)            B A z  =  λ  z .

При этом для всех видов уравнений (1), (2) или (3) искомые значения  λ   всегда являются вещественными.

Вычисляемые матрицы Z, столбцами которых являются собственные векторы z, удовлетворяют соотношениям
              ZT A Z   =  Λ  -  для вещественных матриц
     или   
              ZH A Z  =  Λ  -  для эрмитовых матриц

     при решении уравнений вида (1) или (3), а также соотношениям
Z -1 A Z -T  =  I -  для вещественных матриц
      или   
Z -1 A Z -H  =  I  -  для эрмитовых матриц

      при решении уравнения вида (2).
Здесь Λ - диагональная матрица с расположенными на диагонали собственными значениями, символом T обозначена операция транспонирования, символом H - операция комплексного сопряжения.
Кроме того, матрицы  Z  удовлетворяют соотношениям

ZT B Z   =   I  -  для вещественных матриц

      или   
                 ZH B Z   =   I   -  для эрмитовых матриц,

при решении уравнений (1) или (2), а также соотношениям

ZT B -1 Z   =   I  -  для вещественных матриц
      или   
ZH B -1 Z   =   I  -  для эрмитовых матриц

      при решении уравнения  (3).
Подробнее, вычисление состоит из следующих трех этапов.
1.	Каждый из указанных видов обобщенной проблемы собственных значений (1), (2) или (3) может быть приведен к стандартной (линейной) симметричной проблеме собственных значений посредством использования разложения Холецкого для матрицы B:

B = L LT   или   B = UT U  -  для вещественного симметричного случая,

B = L LH   или   B = UH U  -  для случая эрмитовых матриц.

      Здесь матрицы L и U – нижняя и верхняя треугольные матрицы соответственно. Тогда с учетом разложения B = L LT  получим

A z  =  λ B z     = >     (L -1 A L -T) (LT z)  = λ (LT z)  .

Собственные значения задачи A z = λ B z  являются также собственными значениями задачи  C y = λ y,  где  симметричная матрица  C представляется в виде С= L -1 A L -T ,  а y = LT z .

2.	Выполняется решение полученной стандартной (линейной) проблемы собственных значений, с использованием целевых подпрограмм имеющихся в разделе комплекса для решения линейной симметричной проблемы собственных значений.

3. Восстановление собственных векторов z для обобщенной проблемы собственных значений  из полученных собственных векторов y для линейной проблемы собственных значений  может быть выполнено с использованием соответствующих подпрограмм из пакета BLAS уровня 2 или 3.  

     В таблице ниже показано, как каждая из трех видов проблемы собственных значений может быть сведена к стандартной форме  C y = λ y  и как собственные векторы  z исходной проблемы могут быть получены из собственных векторов  y  стандартной (линейной) проблемы. В таблице представлены только вещественные случаи. Для комплексных случаев транспонированные матрицы должны быть заменены на комплексно-сопряженные транспонированные.

            Вид                       Разложение               Вид                                     Пересчет                 
            проблемы             матрицы B                матрицы C                          векторов z                 
                --------------------------------------------------------------------------------------------------------------------
                  Az = λ  Bz               B = LLT                  C = L -1 AL -T                          z = L –T y
                                B = UT U                 C = U -T AU -1                         z = U -1 y
                 ABz =  z                   B = LLT                  C = LT AL                                z = L -T y
                                B = UT U                 C = UA UT                              z = U -1 y 
                  BAz = λ  z               B = LLT                  C = LT AL                                z = L y
                                B = UT U                 C = UA UT                              z = UT y 

3.	Практические сведения по использованию параллельных программ

3.1.	Предварительные действия, необходимые при обращении к программам

С учетом сказанного в п. 2.2, прежде всего необходимо начинать с определения конкретной решетки процессов. В программах комплекса PARALG используются следующие обозначения:
NPROW - число строк решетки процессов,

NPCOL - число столбцов решетки процессов,

NP = NPROW*NPCOL - общее число процессов решетки.

Для идентификации каждого процесса используется пара (MYROW, MYCOL):
MYROW - номер строки процесса в решетке,

MYCOL - номер столбца процесса в решетке.

      При этом должны быть выполнены неравенства

                   0 ≤  MYROW < NPROW,
                   0 ≤  MYCOL  < NPCOL.

Еще одна процедура, которую необходимо проделать пользователю перед обращением к одной из целевых программ комплекса, это разделить исходную(ые) матрицу(ы) на блоки и распределить эти блоки по определенным правилам по всем процессам решетки.
Это необходимо проделать по той причине, что программы ориентированы на более эффективную блочную реализацию операций с матрицами. Поэтому при обращении к подпрограммам комплекса необходимо задавать, в частности, величину указанных блоков, т.е.:
MB - число строк матрицы в блоке;
NB - число столбцов матрицы в блоке.
Описание конкретных правил распределения матриц по параллельным процессам приводится далее в п. 3.2, а также в п. 2.4.
Проделав, указанные действия, пользователь может вызвать одну из программ комплекса, передав ей через параметры всю информацию об исходных матрицах, а также о том, как распределены их элементы в локальную память каждого из процессов решетки.
После выполнения требуемых вычислений и получения результатов пользователю необходимо освободить использованную решетку процессов и осуществить выход из пакета BLACS.
Ниже подробнее описываются общие обязательные действия программиста при обращении к любой целевой программе Комплекса.
3.1.1.	Инициализация решетки процессов

После присвоения конкретных значений параметрам решетки NPROW и NPCOL инициализация решетки процессов проще всего производится посредством обращения к подпрограмме SL_INIT из разряда служебных ( tool  routines). Она состоит из обращений к необходимым подпрограммам пакета BLACS.  
Эта подпрограмма инициализирует выбранную решетку процессов, упорядочивая процессы по строкам, а также присваивает служебному параметру ICTXT некоторое целое значение. Это значение закрепляет за выбранной пользователем решеткой процессов статус конкретной области действия массовых операций передачи сообщений. Параметр ICTXT называется указателем контекста, системным контекстом, контекстом BLACS'а или просто контекстом.
Пользователь, работая с выбранной решеткой, не должен производить никаких действий с параметром ICTXT. Необходимо только передавать его в качестве входного параметра при обращении к другим подпрограммам пакета BLACS, а также к подпрограммам комплекса, куда передаются сведения об обрабатываемых матрицах.
Обращение к подпрограмме SL_INIT выглядит следующим образом.
          CALL SL_INIT ( ICTXT, NPROW, NPCOL) ,
     где
ICTXT     - специфицирует контекст BLACS'а (глобальный выходной параметр, тип:
                   целый);
NPROW  - число строк, в создаваемой решетке процессов (глобальный входной
                   параметр, тип: целый);
NPCOL   - число столбцов, в создаваемой решетке процессов (глобальный входной 
                   параметр, тип: целый).
Программы комплекса написаны с расчетом на использование в стиле SPMD (Single Program Multiple Data). Это означает, что одна и та же программа пользователя будет загружена и одновременно выполняться всеми процессами решетки. Поэтому при необходимости выполнения некоторых действий только на некоторых из этих процессов требуется определить координаты процесса, на котором выполняется данный экземпляр программы. Как уже указывалось, они обычно обозначаются идентификаторами MYROW и MYCOL. Например, может быть необходимо произвести распечатку входных параметров, одинаковых для всех процессов. Чаще всего это делают с процесса с координатами (0, 0).
Для определения координат процесса используют подпрограмму BLACS_GRIDINFO пакета BLACS, обращение к которой имеет следующий вид. 
          CALL BLACS_GRIDINFO ( ICTXT, NPROW, NPCOL, MYROW, MYCOL) ,
      где MYROW - номер строки данного процесса в решетке процессов,
                                0 ≤  MYROW < NPROW (глобальный выходной параметр);

              MYCOL  - номер столбца данного процесса в решетке процессов,
                                 0 ≤ MYCOL < NPCOL (глобальный выходной параметр).

      Три первых параметра описаны выше, только в данном случае параметр ICTXT - входной.

3.1.2. Распределение матриц по решетке процессов

Поскольку каждый процесс производит операции над своей частью исходной глобальной матрицы, до обращения к программам комплекса необходимо распределить по определённым правилам матрицы A и B по решетке процессов. Это является обязанностью пользователя (подробности этого процесса описаны в п. 2.4, а также в п. 3.2).
При этом каждой глобальной матрице, распределенной по решетке процессов, должен быть приписан дескриптор, содержащий информацию о том, как именно было произведено это распределение. Дескриптор представляет собой одномерный массив из 9 или 7 элементов целого типа. Подробное описание дескрипторов можно найти в п. 2.3.2. Здесь мы покажем, как инициализировать дескриптор DESCA глобальной матрицы A. 
Если матрица A - плотная, то проще всего инициализировать дескриптор посредством обращения к служебной подпрограмме DESCINIT.
Ниже приводится обращение к этой подпрограмме при инициализации дескриптора DESCA. Подобные обращения можно найти, например, в текстах тестовых примеров к программам TDSYGV11.f и TZHEGV11.f вычисления всех собственных значений вещественных симметричных или эрмитовых матриц в обобщенной проблеме собственных значений. 
          CALL DESCINIT(DESCA, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA, INFO)
Этот вызов подпрограммы DESCINIT эквивалентен следующим операторам присваивания:
DESCA(1) = 1
DESCA(2) = ICTXT
DESCA(3) = M
DESCA(4) = N
DESCA(5) = MB
DESCA(6) = NB
DESCA(7) = RSRC
DESCA(8) = CSRC
DESCA(9) = MXLLDA

Здесь первый элемент обозначает тип дескриптора для плотной матрицы, второй - обсуждавшийся выше служебный параметр подпрограмм BLACS'а (контекст).
Другие параметры имеют следующий смысл:
M, N          - число строк и столбцов глобальной матрицы  A ( M , N );
MB, NB     - число строк и столбцов в блоках, на которые разбивается матрица A;
RSRC         - координаты процесса в решетке, куда размещается пользователем первый
CSRC           элемент распределяемой глобальной матрицы;
         MXLLDA - ведущая локальная размерность матрицы A.

         INFO        - целая переменная, служащая для сообщения о результате работы программы 
                             (выходной параметр).

Простой пример отображения глобальной матрицы на решетку процессов можно найти, в частности, в тестовом примере к подпрограмме вычисления собственных значений вещественных симметричных матриц PDSYGV1.f (см. TDSYGV1.f) в обобщенной проблеме собственных значений. В нем эти действия выполняются в подпрограмме с именем MATINIT (или PDMATINIT), которая является частью тестового примера.

3.1.3	Освобождение решетки процессов

После того как необходимые вычисления на решетке процессов были завершены, требуется освободить решетку процессов посредством вызова подпрограммы BLACS_GRIDEXIT. Когда завершены все вычисления, выход из пакета BLACS должен быть осуществлен посредством обращения к подпрограмме BLACS_EXIT.
Типичный пример фрагмента программы, включающий эти действия: 
    CALL BLACS_GRIDEXIT(ICTXT)
CALL BLACS_EXIT(0)

3.2.	Программирование распределения матриц по решетке процессов

В настоящем учебном пособии рассматриваются три способа реализации на языке ФОРТРАН рассылки элементов исходных глобальных матриц по параллельным процессам, поскольку каждый из них удобнее использовать в трех разных ситуациях.
1.	Первый способ удобен только для небольших матриц и требует от пользователя самостоятельного размещения элементов исходных матриц в локальной памяти параллельных процессов в соответствии с общими правилами, описанными в п. 2.4.  При этом используются обычные операторы присваивания языка ФОРТРАН. Такой способ достаточно трудоемок и может использоваться для целей обучения и усвоения общих правил распределения матриц (см. п. 2.4).

2.	Второй способ удобнее как для небольших матриц, так и для матриц среднего размера с регулярной структурой или для матриц любого размера, элементы которых задаются математическими формулами.

Этот способ гораздо проще для пользователя, поскольку не требует самостоятельного определения местоположения элементов матриц в локальной памяти каждого из параллельных процессов. Он опирается на использование служебной подпрограммы PDELSET, которая сама (в соответствии с упомянутыми общими правилами, см. п.2.4.1) заносит элемент с индексами  I, J исходной матрицы  A ( I, J ) в определенное место локальной памяти того или другого из параллельных процессов.
На ФОРТРАНе эта процедура сводится к написанию некоторого числа циклов, внутри которых стоит обращение к подпрограмме PDELSET.
Фортранная подпрограмма с реализацией такого способа распределения приводится в конце Приложения 1.
Подробный разбор приводимого в Приложениях примера дается в п. 3.4.
3.	Третий способ может быть использован для обработки больших исходных матриц, численное значение элементов которых хранится в виде файлов на внешней памяти (наиболее частый на практике случай).

Для реализации третьего способа можно воспользоваться специальной служебной подпрограммой PDLAREAD, последовательно считывающей из внешнего файла элементы исходной матрицы и распределяющей их по установленным правилам в локальную память того или другого из параллельных процессов.
Однако пользователь может написать свою аналогичную подпрограмму, если какие-либо детали размещения исходной матрицы в файле отличаются от приводимого примера.
Кроме того, возможно обращение к служебной подпрограмме PDLAWRITE для записи полученных результатов в файл. Тексты этих подпрограмм можно найти, перейдя по гиперссылкам в конце Интернет-страницы с адресом 
http://num-anal.srcc.msu.ru/par_prog/instr.htm.

3.3.	O документировании и примерах использования параллельных программ

	Правила наименований подпрограмм

Каждая подпрограмма комплекса PARALG имеет имя, начинающееся с буквы P. Допущено нарушение стандарта языка Фортран-77, так как имена подпрограмм могут быть длиной более шести символов; кроме того, в именах служебных подпрограмм (tool routines) допускается использование символа подчеркивания “ _”.
Все целевые и базовые подпрограммы имеют имена в виде последовательности шести или семи символов: PXYYZZ или PXYYZZZ, где второй символ X может принимать следующие значения, указывающие на основной тип обрабатываемых данных:
S    REAL
D    DOUBLE PRECISION
C    COMPLEX
Z    DOUBLE COMPLEX 

   Следующие две буквы YY указывают на то, какого вида матрица обрабатывается подпрограммой (или вид самой главной из участвующих матриц). Большинство из этих двухсимвольных сочетаний относятся одновременно как к вещественным, так и к комплексным матрицам, некоторые же относятся конкретно к одному из этих типов матриц. Приняты следующие мнемонические правила для указания вида обрабатываемых матриц:

DB (general band)          - ленточная матрица общего вида, для которой не 
                                         требуется выбор ведущего элемента
DT (general tridiagonal) - трехдиагональная матрица общего вида, для которой не
                                          требуется выбор ведущего элемента
GB (general band)          - ленточная матрица общего вида
          GE (general)                   - матрица общего вида (т.е. не симметричная, иногда 
                                                    прямоугольная) 
           GG (general matrices, generalized problem) - матрица общего вида, обобщенная
                                                                                     проблема собственных значений (т.е.
                                                                                     пара матриц общего вида)
HE (complex Hermitian) - комплексная эрмитова матрица
OR (real orthogonal)       - вещественная ортогональная матрица
PB (symmetric or Hermitian positive definite band) - симметричная или эрмитова 
                                                                                      ленточная положительно 
                                                                                      определенная
PO (symmetric or Hermitian positive definite) - симметричная или эрмитова
                                                                              положительно определенная
PT (symmetric or Hermitian positive definite tridiagonal) – симметричная или 
                                                                                            эрмитова трехдиагональная 
                                                                                            положительно определенная 
ST (real symmetric tridiagonal) - вещественная симметричная трехдиагональная
SY (symmetric)                         - симметричная матрица
TR (triangular or in some cases quasi-triangular) – трехдиагональная матрица или в 
                                                               некоторых случаях почти трехдиагональная
TZ (trapezoidal) – трапециевидная матрица
UN (complex unitary) - комплексная унитарная матрица

Последние два символа ZZ или три ZZZ символа в имени подпрограммы указывают ее назначение. Например, окончание TRF означает подпрограммы, которые выполняют факторизацию матриц.
Имена вспомогательных подпрограмм подчиняются тем же правилам, за исключением того, что третьим и четвертым символом YY обычно являются символы LA (например, PDLASCL или PZLARFG).
	Описания программ и примеры их использования

        Документация параллельных программ по линейной алгебре включает в себя:
- информацию общего характера (об организации и структуре Комплекса и об общих правилах использования программ и размещения данных); 
- систематический каталог целевых программ; 
- описания целевых программ (т.е. инструкции по их использованию).  
Документация общего характера приводится как в рамках настоящего учебного пособия, так и в соответствующих подразделах раздела «Комплекс программ по линейной алгебре для вычислений на распределенной памяти» Научно - образовательного Интернет -ресурса НИВЦ МГУ по численному анализу» (http://num-anal.srcc.msu.ru/). Там же приводится систематический каталог целевых программ в виде таблиц, содержащих, в частности, гипер-ссылки на описания программ.
Описания программ составлялись во многом по аналогии с тем, как это было ранее разработано и реализовано в Библиотеке фортранных программ по численному анализу НИВЦ МГУ [3], которая состоит из  последовательных алгоритмов.  
Такая четкая форма документации, состоящая из нескольких обязательных смысловых разделов, заслужила широкое одобрение среди многочисленных пользователей Библиотеки численного анализа, в течение многих лет эксплуатировавших ее в разных организациях и на разных вычислительных системах.  
Одним из самых важных разделов в описании каждой целевой программы является его последний раздел под названием «Пример использования». Содержание этого раздела на примере матриц небольшого размера наглядно демонстрирует пользователю, как будут распределены их элементы по параллельным процессам и как такое распределение можно реализовать средствами языка ФОРТРАН и служебных подпрограмм. А именно, приводится фрагмент головной фортранной программы (теста), с подпрограммой рассылки элементов конкретной матрицы по параллельным процессам одним из трех способов, рассматриваемых в настоящем пособии (см. п.3.2.). Указанные подпрограммы рассылки имеют в примерах либо имя MATINIT, либо PDMATINIT.
Обозначения в примерах по использованию программ Комплекса

В тестовых примерах и примерах по использованию программ комплекса приняты следующие единообразные обозначения величин, задающих параметры исходной задачи и выбранный способ распараллеливания (указаны в алфавитном порядке).
  Имя                                            Определение

CSRC          - номер столбца решетки процессов, в который распределяется первый
                     столбец исходной глобальной матрицы

DESCA       - дескриптор исходной глобальной матрицы  A

DESCB       - дескриптор исходной глобальной матрицы (вектора) B (или правых
                     частей при решении систем линейных уравнений, или при решении 
                     обобщённой проблемы собственных значений)

ICTXT       - контекст BLACS'а, связанный с выбранной решеткой процессов
                    (внутренний параметр комплекса; служит для обеспечения согласованной
                    работы подпрограмм комплекса и устанавливается пользователем 
                    обращением к соответствующей подпрограмме BLACS'а в начале решения
                    задачи, см. пп.3.1. и 2.3.2.)

M              - число строк исходной глобальной матрицы  A

MB           - число строк в блоке, на которые делится исходная матрица  A

MXLLDA - максимальная ведущая локальная размерность локальной части
                    матрицы  A среди всех локальных частей на всех процессах решетки

MXLLDB - максимальная ведущая локальная размерность локальной части
                    матрицы (вектора)  B среди всех локальных частей на всех процессах 
                    решетки

MXLOCC - максимальное число столбцов в локальных частях матрицы  A,
                     которые принадлежат процессам в столбцах решетки процессов

MXLOCR - максимальное число строк в локальных частях матрицы  A, которые
                    принадлежат процессам в строках решетки процессов

MXRHSC - максимальное число столбцов в локальных частях матрицы (вектора)
                    правых частей  B, которые принадлежат процессам в столбцах решетки
                    процессов (если  B - вектор, то MXRHSC = 1)
MYCOL   - координата (номер) столбца вызываемого процесса в решетке процессов

MYROW - координата (номер) строки вызываемого процесса в решетке процессов

N             - число столбцов исходной глобальной матрицы  A (и число строк
                  глобальной матрицы (вектора) правых частей  B при решении систем
                 уравнений)

NB         - число столбцов в блоке, на которые делится исходная матрица  A

NBRHS - число столбцов в блоке, на которые делится глобальная матрица (вектор)
                правых частей системы уравнений  B (если  B - вектор, то NBRHS =1)

NPCOL - число столбцов в решетке процессов

NPROW - число строк в решетке процессов

NRHS    - число столбцов исходной глобальной матрицы (вектора) правых частей
                 системы уравнений  B (если  B - вектор, то NBRHS = 1)

RSRC    - номер строки решетки процессов, в которую распределяется первая
                строка исходной глобальной матрицы

Упрощения, принятые в примерах по использованию программ комплекса

Опишем далее некоторые упрощения и ограничения, которые были сделаны при подготовке примеров использования подпрограмм комплекса, для того чтобы упростить пользователю освоение правил работы с этими подпрограммами.
1.	В примерах предполагается, что RSRC = CSRC = 0; это означает, что исходные матрицы  A (и  B) распределяются по решетке процессов начиная с процесса (0, 0). На самом деле, любой из процессов решетки может быть выбран в качестве начального при распределении матриц.

2.	В примерах предполагается, что задача решается при использовании всей исходной распределенной матрицы  A (хотя подпрограммы допускают, чтобы задачи решались с матрицей, являющейся подматрицей sub (A) исходной матрицы  A). При указанном предположении параметры, задающие первый элемент подматрицы  A, полагаются равными 1, т.е.  IA = JA = IB = JB = 1.

3.	При распределении исходных матриц  A по решетке процессов может использоваться подпрограмма MATINIT, которая заносит с помощью операторов присваивания в соответствующие локальные массивы на каждом из процессов решетки распределенные на него элементы исходных матриц. Упрощает эту процедуру использование служебной подпрограммы PDELSET в подпрограмме PDMATINIT (см. п.  3.2. и конец Приложения 1).

4.	В примерах предполагается, что ведущая локальная размерность локальных частей матрицы  A является одинаковой для всех процессов, принадлежащих к одним и тем же строкам решетки процессов (т.е. имеющих одинаковые координаты строк решетки).  Переменная MXLLDA равна максимальной ведущей локальной размерности для матрицы A (обозначаемой LLD_A) из всех строк решетки процессов. Однако в общем случае локальные ведущие размерности локальных частей матриц могут отличаться друг от друга для всех процессов решетки.

5.	В примерах часто рассматриваются исходные матрицы  небольших размеров. При этом размеры блоков, на которые разбиваются эти матрицы, выбираются небольшими.  Например,  MB = NB = 2. Однако это не является самым лучшим размером блоков для практических задач. Для достижения лучшей производительности более подходящим размером блоков будет  MB = NB = 32 или  MB = NB = 64.

3.4. Пример использования программы решения обобщенной проблемы собственных значений для симметричных матриц

     Рассмотрим подробнее пример по использованию подпрограммы PDSYGV2, которая предназначена для вычисления всех собственных значений и собственных векторов в обобщенной проблеме собственных значений для вещественных симметричных матриц общего вида  A и  B. Данный пример приводится в описании этой подпрограммы в разделе «Пример использования» (см. Приложение 1).  

     Пример демонстрирует решение обобщенной проблемы собственных значений для симметричных матриц  A и  B размера 5 x 5. Ниже приводится их изображение в символическом и числовом виде.

a11   a12   a13   a14   a15        10    2     3    1    1
a21   a22   a23   a24   a25         0   12    1     2    1
a31   a32   a33   a34   a35   =    0    0    11    1   -1
a41   a42   a43   a44   a45         0    0     0     9    1
a51   a52   a53   a54   a55         0    0     0     0  15

b11   b12   b13   b14   b15        12    1   -1     2    1
b21   b22   b23   b24   b25         0   14    1    -1    1
b31   b32   b33   b34   b35   =    0    0    16   -1    1
b41   b42   b43   b44   b45         0    0     0    12   -1
b51   b52   b53   b54   b55         0    0     0     0   11

Порядок матриц N = 5.

     Сначала мы должны задать решетку процессов и распределить по ней исходные глобальные матрицы  A и  B.

     Как говорилось в пп.3.1, 2.4.1, базовые подпрограммы комплекса реализуют специальные блочные алгоритмы обработки матриц, позволяющие повысить скорость решения задачи. При этом распределение плотных матриц по параллельных процессам происходит блочно-циклическим образом (см. пп.2.4.1, 2.4.2 и п.3.2 ).

     Пусть матрицы разбиваются на блоки размером  MB x NB, где  MB = NB = 2, следующим образом.

10    2   |   3    1   |   1
           0    12  |   1    2   |   1
           ----------------------------
           0    0    |   11  1   |  -1
           0    0    |   0    9   |   1
          ---------------------------
           0    0    |   0    0   |  15

12    1   | -1    2   |   1
           0    14  |  1     -1  |   1
           ----------------------------
           0    0    |   16  -1  |   1
           0    0    |   0    12 |  -1
          ---------------------------
           0    0    |   0    0   |   11

     Пусть пропуск программы осуществляется на 4 процессах, которые образуют решетку 2 на 2, т.е. NPROW = 2 и NPCOL = 2 (см. пп.2.2, 3.1 ).

     После распределения полученных блоков по процессам решетки (по принятым правилам, которые описаны в пп.2.4.1, 2.4.2), имеем следующую картинку распределения элементов матрицы  A по процессам (в символическом и числовом виде). Для матрицы B распределение выполняется аналогично. Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.

                           0                         1                         
---------------------------------------------------
           a11   a12      a15         a13   a14
           a21   a22      a25         a23   a24
  0       -------------------------------------------
           a51   a52      a55         a53   a54
----------------------------------------------------
           a31   a32       a35         a33   a34
  1       a41   a42       a45         a43   a44

                            0                       1                         
----------------------------------------------------
           10     2           1             3        1
             0    12          1             2        1
  0        --------------------------------------------
             0     0           15           0        0
-----------------------------------------------------
             0     0           -1           11       1
  1         0     0            1             0        9

Из последнего рисунка видно, что процесс с координатами (0, 0) содержит локальный массив (часть матрицы  A) размера (3, 3), процесс (0, 1) - локальный массив размера (2, 3), процесс (1, 0) - локальный массив размера (3, 2), процесс (1, 1) - локальный массив размера (2, 2).
Для небольших матриц возможен способ инициализации элементов распределенных матриц с помощью обычных операторов присваивания, каждый из которых выполняется только на одном из параллельных процессов. Такую подпрограмму пользователь может написать только заранее представив себе распределение элементов матриц в соответствии с изображенными рисунками, а также учитывая, что в языке Фортран двумерные массивы хранятся в памяти по столбцам.
Такой способ инициализации распределенных матриц явно не годится для матриц больших размеров. В этой связи в Приложении 1 приводится более простой способ инициализации матриц с помощью обращения к служебной подпрограмме PDELSET, которая сама определяет (в соответствии с принятыми правилами блочно-циклического распределения), на какой из параллельных процессов и в какое место локальной памяти (для локальной части матрицы) следует поставить конкретный элемент глобальной матрицы с глобальными индексами   I, J,  т.е.  A ( I, J ) или B ( I, J ) .
Еще один способ распределения исходной матрицы по параллельным процессам (возможно, имеющий наибольшее практическое применение) – это случай считывания элементов исходной матрицы из внешнего файла. Такое считывание и последующее распределение ее элементов по решетке параллельных процессов выполняется с помощью служебной подпрограммы PDLAREAD. Такой способ распределения матрицы является наиболее автоматизированным. Он требует от пользователя только правильного выделения локальной памяти под локальные части исходных матриц, массива собственных векторов и рабочего массива.
Вычислить размеры локальной части распределенного двумерного массива можно с помощью служебной подпрограммы NUMROC, описание функций которой приводится в п.2.3.2.
3.5.	Общий список программ комплекса для решения обобщенной алгебраической проблемы собственных значений для симметричных и эрмитовых матриц.

В п. 2.5 настоящего пособия были перечислены алгоритмы, используемые в комплексе при решении обобщенной алгебраической проблемы собственных значений для симметричных и эрмитовых матриц.

Приведем здесь список всех целевых программ, включенных в настоящее время в раздел комплекса «Решение обобщенной проблемы собственных значений для симметричных и эрмитовых матриц».
         PDSYGV1 - вычисление всех собственных значений в обобщенной проблеме 
         PZHEGV1  собственных значений для вещественных симметричных (эрмитовых)
                             матриц A и B

PDSYGV2 - вычисление всех собственных значений и собственных векторов в
PZHEGV2   обобщенной проблеме собственных значений для вещественных 
                     симметричных (эрмитовых) матриц A и B

PDSYGV3 - вычисление собственных значений, принадлежащих заданному
PZHEGV3   интервалу (VL, VU), в обобщенной проблеме собственных значений 
                     для  вещественных симметричных (эрмитовых) матриц A и B

PDSYGV4 - вычисление собственных значений, принадлежащих заданному
PZHEGV4   интервалу (VL, VU), и собственных векторов, соответствующих этим
                     собственным значениям, в обобщенной проблеме собственных значений
                    для вещественных симметричных (эрмитовых) матриц A и B

PDSYGV5 - вычисление собственных значений, принадлежащих заданному
PZHEGV5   интервалу индексов (IL, IU), в обобщенной проблеме собственных 
                     значений для вещественных симметричных (эрмитовых) матриц A и B

PDSYGV6 - вычисление собственных значений, принадлежащих заданному
PZHEGV6   интервалу индексов (IL, IU), и собственных векторов, соответствующих 
                    этим  собственным значениям, в обобщенной проблеме собственных 
                    значений для вещественных симметричных (эрмитовых) матриц A и B

     Дадим здесь некоторые пояснения к реализации указанных в п.2.5. алгоритмов базовыми подпрограммами комплекса.

     Разложение Холецкого матрицы B выполняется с помощью базовой подпрограммы PDPOTRF. Используя это разложение, приведение обобщенной проблемы собственных значений к стандартной линейной форме выполняется с помощью базовых подпрограмм PDSYNGST (PDHEGST).

     Далее используется, необходимая для конкретного случая, одна из целевых подпрограмм комплекса для решения линейной проблемы собственных значений: PDSYEV1, PZHEEV1, PDSYEV2, PZHEEV2, PDSYEV3, PZHEEV3, PDSYEV4, PZHEEV4, PDSYEV5, PZHEEV5, PDSYEV6, PZHEEV6.

При необходимости вычисления собственных векторов, восстановление собственных векторов для обобщенной проблемы собственных значений по собственным векторам линейной проблемы собственных значений производится с помощью базовых подпрограмм PDTRSM (для уравнений вида (1) или (2)) или PDTRMM (для уравнений вида (3)).

3.6.	 Запуск программ комплекса в ОС Linux на суперкомпьютере “Чебышев”

Комплекс программ PARLAG доступен пользователям суперкомпьютера СКИФ “Чебышев” в НИВЦ МГУ в виде откомпилированных библиотек. Пакеты (BLAS, BLACS, ScaLAPACK и др.), модули которых используются программами комплекса, входят в состав MKL-библиотеки (Intel Math Kernel Library), установленной на этом суперкомпьютере, и должны быть подсоединены при получении исполнимой программы.  Объектные модули самого комплекса берутся из библиотеки libparalg.a, указанной в последней строке приводимого ниже скрипта bld-tcl, который предназначен для компиляции и получения исполнимой программы для ее последующего запуска на счет в параллельном режиме на суперкомпьютере “Чебышев”.

#! /bin/sh
EXEDIR=.
EXE=”<имя исполнимой программы>”
OBJ=”<список имен объектных модулей>”
F77OPTS=”-O2”
F77=mpif77
#
$F77 -c $F77OPTS <имя исходного фортранного модуля>
#
echo Linking executable $EXE
#
$F77 -o $EXEDIR/$EXE $OBJ \
            -Wl,--start-group -lmkl_scalapack_lp64 -lmkl_blacs_lp64 \
            -lmkl_intel_lp64 \
            -lmkl_sequential -lmkl_core \
            -Wl,--end-group \
            -lpthread -libparalg.a

Здесь:
•	<имя исполнимой программы> - указываемое пользователем имя, которое будет присвоено полученному исполнимому модулю (например, tdsygv21.e);
•	<список имен объектных модулей> - список имен всех объектных модулей, которые не вызываются из библиотек (например, tdsygv21.o);
•	$F77 -c $F77OPTS <имя исходного фортранного модуля> - команда компиляции (которых может быть несколько) фортранного модуля (например, $F77 -c $F77OPTS tdsygv21.f).

В скрипте bld-tcl подключаются все пакеты из MKL-библиотеки, содержащие все вызываемые подпрограммы, используемые при работе параллельных подпрограмм комплекса. Полученная исполнимая программа может быть затем запущена на счет с помощью команды mpirun с параметрами.
Например, если требуется запустить программу tdsygv21.e на шести процессорах, то команда запуска имеет вид

mpirun -np 6 -maxtime 1 ./tdsygv21.e

В тестах к программам выдача исходных матриц и полученных результатов осуществляется с помощью обращения к подпрограмме PDLAPRNT из пакета ScaLAPACK(TOOLS). Эта подпрограмма распечатывает всю распределенную матрицу (вектор) полностью по столбцам по одному элементу в строке, получая элементы со всех участвующих в вычислениях процессоров.
Обращение к программам комплекса может осуществляться с помощью подпрограмм автоматизации доступа к ним. Описание этих подпрограмм доступно по адресу 
http://num-anal.srcc.msu.ru/par_prog/org/avt.htm.
Тогда конкретный вид скрипта может быть проиллюстрирован следующим примером для запуска на счет программы PDSYGV2:
#! /bin/sh
EXEDIR=.
EXE=”tcl_sygvz.e”
OBJ=”tcl_sygvz.o cal_sygvz.o pdsygv2.o pdlaprnt.o pdlaread.o pdlawrite.o”
# 
F77OPTS=”-O2” 
F77=mpif77
#
         $F77 -c $F77OPTS tcl_sygvz.f
#
         $F77 -o $EXEDIR/$EXE $OBJ \
           -Wl,--start-group -lmkl_scalapack_lp64 -lmkl_blacs_lp64 \
           -lmkl_intel_lp64 \
           -lmkl_sequential -lmkl_core \
           -Wl,--end-group \
           -lpthread -libparalg.a

Чтобы не указывать все объектные модули подпрограмм комплекса PARALG, которые необходимо вызвать для решения конкретной задачи, достаточно указать в последней строке скрипта объектную библиотеку libparalg.a. Тогда строка с объектными модулями будет содержать только объектный модуль головной программы OBJ=“tcl_sygvz.o”.
Здесь подпрограммы pdlaread.o и pdlawrite.o - входящие в состав комплекса служебные подпрограммы для считывания из файла исходных матриц и распределения их по параллельным процессам и записи в файл полученного вектора решения.

3.7.	 Запуск программ в ОС Linux на суперкомпьютере “Ломоносов-1”

Комплекс программ PARLAG доступен пользователям суперкомпьютера “Ломоносов-1” в НИВЦ МГУ в виде откомпилированных библиотек. Пакеты (BLAS, BLACS, ScaLAPACK и др.), модули которых используются программами комплекса, входят в состав MKL-библиотеки (Intel Math Kernel Library ), установленной на этом суперкомпьютере(версия 11.2.0), и должны быть подсоединены при получении исполнимой программы. Объектные модули самого комплекса берутся из библиотеки libparalg.a, указанной в последней строке приводимого ниже скрипта bld, который предназначен для компиляции и получения исполнимой программы для ее последующего запуска на счет в параллельном режиме на суперкомпьютере “Ломоносов-1”.
#! /bin/sh
EXEDIR=.
EXE=”<имя исполнимой программы>”
OBJ=”<список имен объектных модулей>”
F77OPTS=”-O2”
F77=mpif77
         #
$F77 -c $F77OPTS <имя исходного фортранного модуля>
        #
echo Linking executable $EXE
        #
$F77 -o $EXEDIR/$EXE $OBJ \
                             -lmkl_scalapack_lp64 -lmkl_intel_lp64 \
                             -lmkl_sequential -lmkl_core \
                             -lmkl_blacs_intelmpi_lp64  \
                             -lpthread -lm -ldl libparalg.a

Здесь:
•	<имя исполнимой программы> - указываемое пользователем имя, которое будет присвоено полученному исполнимому модулю (например, tdsygv21.e);
•	<список имен объектных модулей> - список имен всех объектных модулей, которые не вызываются из библиотек (например, tdsygv21.o);
•	$F77 -c $F77OPTS <имя исходного фортранного модуля> - команда компиляции (которых может быть несколько) фортранного модуля 
      (например,  $F77 -c $F77OPTS tdsygv21.f).

По правилам работы на Ломоносов-1 перед получением исполнимой программы (т.е.  перед запуском скрипта bld) следует установить необходимое окружение. Для этого необходимо выполнить команду 
module  add  intel/15.0.090  impi/5.0.1-ofa  slurm/2.5.6  mkl/11.2.0,
которая загрузит необходимые версии компилятора с Фортрана, библиотеки MPI и MKL-библиотеки.
В скрипте bld подключаются все пакеты из MKL-библиотеки, содержащие все вызываемые подпрограммы, используемые при работе параллельных программ комплекса. Перед запуском на счет (постановкой в очередь) необходимо перейти в подкаталог _scratch и перенести туда библиотеку комплекса libparalg.a.  Полученная исполнимая программа может быть затем запущена на счет с помощью команды sbatch с параметрами.
Например, если требуется запустить программу tdsygv21.e на восьми процессорах, то команда запуска может иметь вид 
sbatch  -n 8  -p test  -t 1  -o tdsygv21.res  impi  ./tdsygv21.e 
Здесь указано, что программа должна выполняться в тестовой очереди, заказанное для выполнения время равно одной минуте, а результаты следует поставить в файл tdsygv21.res.
В тестах к программам выдача исходных матриц и полученных результатов осуществляется с помощью обращения к подпрограмме PDLAPRNT из пакета ScaLAPACK(TOOLS). Эта подпрограмма распечатывает всю распределенную матрицу (вектор) полностью по столбцам по одному элементу в строке, получая элементы со всех участвующих в вычислениях процессоров.
Чтобы не указывать все объектные модули подпрограмм комплекса PARALG, которые необходимо вызвать для решения конкретной задачи, достаточно указать в последней строке скрипта объектную библиотеку libparalg.a. Тогда строка с объектными модулями будет содержать только объектный модуль головной программы OBJ=“tdsygv21.o”.

Приложение 1

Пример описания целевой подпрограммы для вычисления всех собственных значений и собственных векторов вещественных симметричных матриц
 в обобщенной проблеме собственных значений

Подпрограмма:  PDSYGV2

Назначение

Вычисление всех собственных значений и собственных векторов в обобщенной проблеме собственных значений для вещественных симметричных матриц A и B 
Математическое описание 
Данная подпрограмма реализует алгоритм вычисления всех собственных значений и собственных векторов для задач вида
                              (1)                     Аx  = λ Bx,

                              (2)                     АBx  = λ x,

(3)	                   BАx  = λ x,

где  A  и  B - вещественные симметричные матрицы и матрица B положительно определена.

При помощи разложения Холецкого для матрицы B вида В = LLT, где L – нижняя треугольная матрица, исходные задачи (1), (2) или (3) приводятся к стандартной задаче  вида Cy = λ y, где

для (1):  C  =  L-1AL-T,  а   y  =  L-Tx .

для (2):  C  =  LTAL   ,    а   y  =  L-Tx .

для (3):  C  =  LTAL   ,    а   y  =  Lx .

Стандартная задача решается путем приведения симметричной матрицы C к трехдиагональному виду и вычисления ее собственных значений λ и собственных векторов y неявным  QL - или QR - алгоритмом.

Восстановление собственных вектоpов  x  исходных задач (1), (2) или (3) из соответствующих вектоpов  у  стандартной задачи осуществляется по формулам

x   =   L-T y      для   (1) или (2)  ,

x   =   L y         для   (3) .

Литература:
http://num_anal.srcc.msu.ru/par_prog/
Дж.Х. Уилкинсон, Алгебраическая проблема собственных значений. M.: Наука, 1970.

Использование
CALL  PDSYGV2 ( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, W, Z, IZ,
                                  JZ, DESCZ, WORK, LWORK, IFAIL, INFO)

Параметры 

IBTYPE - определяет вид обобщенной проблемы собственных значений:
                 если IBTYPE = 1 -  вид (1);
                 если IBTYPE = 2 -  вид (2);
                 если IBTYPE = 3 -  вид (3);
                 (глобальный входной параметр, тип целый);

UPLO     - если UPLO = ' U ' - задается верхний треугольник матрицы A; 
                  если UPLO = ' L ' - задается нижний треугольник матрицы A 
                  (глобальный входной параметр, тип символьный);

N            - порядок распределенной подматрицы sub (A), N > 0 (глобальный входной
                 параметр, тип целый);

A           - массив двойной точности, распределенный по процессам блочно-циклическим
                образом (см. п. 2.4), глобальная размерность которого (N, N), а локальная 
                размерность ( LLD_A, LOCc ( JA + N - 1)); на входе - это локальная часть
                распределенной симметричной подматрицы sub(A); если UPLO = ' U ', то
                используется только верхняя треугольная часть подматрицы sub(A);
                если UPLO = ' L ', то используется только нижняя треугольная часть подматрицы
                sub(A); на выходе - нижний треугольник (при UPLO = ' L ' ) или верхний 
                треугольник (при UPLO = ' U ' ) матрицы sub(A), включая диагональ, не 
                сохраняется (локальный входной параметр и локальный выходной параметр);

IA          - глобальный номер строки матрицы A, который указывает на начало
                подматрицы, которая должна быть обработана (глобальный входной параметр,
                тип целый);

JA          - глобальный номер столбца матрицы A, который указывает на начало
                подматрицы, которая должна быть обработана (глобальный входной параметр,
                тип целый);

DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN);
                     если тип дескриптора одномерный ( DTYPE_A = 501), то DLEN ≥  7; если тип 
                     дескриптора двумерный ( DTYPE_A = 1), то DLEN ≥ 9; дескриптор DESCA
                     содержит информацию о размещении матрицы A в памяти; полное описание
                     DESCA см. в п. 2.3.2 (глобальный и локальный входной параметр, тип целый);

B          - массив двойной точности, распределенный по процессам блочно-циклическим
               образом ( см. п. 2.4 ), глобальная размерность которого (N, N), а локальная
               размерность ( LLD_B, LOCc ( JB + N - 1)); на входе - это локальная часть
               распределенной симметричной подматрицы sub(B);
               если UPLO = ' U ', то используется только верхняя треугольная часть подматрицы
               sub(B); если UPLO = ' L ', то используется только нижняя треугольная часть
                    подматрицы sub(B); на выходе, если INFO ≤ N, то часть матрицы sub(B) 
                    содержит матрицу, замененную на треугольный множитель L из разложения
                     Холецкого  sub(B) = L*LT (локальный входной параметр и локальный выходной
                    параметр);

IB          - глобальный номер строки матрицы B, который указывает на начало
                подматрицы, которая должна быть обработана (глобальный входной параметр,
                тип целый);

JB         - глобальный номер столбца матрицы B, который указывает на начало
               подматрицы, которая должна быть обработана (глобальный входной параметр,
               тип целый);

DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN);
                     если тип дескриптора одномерный ( DTYPE_B = 501), то DLEN ≥ 7; если тип
                     дескриптора двумерный ( DTYPE_B = 1), то DLEN ≥ 9; дескриптор DESCB
                     содержит информацию о размещении B в памяти; полное описание DESCB
                     см. в п. 2.3.2 (глобальный и локальный входной параметр, тип целый);

W          - одномерный массив двойной точности длины N; если  INFO = 0, то собственные 
                значения располагаются в нем по возрастанию (глобальный выходной 
                параметр);

Z           - массив двойной точности с глобальной размерностью ( N, N ) и с локальной
                размерностью ( LLD_Z, LOCc ( JZ + N - 1)); содержит ортонормированные
               собственные векторы симметричной матрицы A (локальный выходной
               параметр);

IZ         - глобальный номер строки матрицы Z, который указывает на начало подматрицы 
               sub(Z) (глобальный входной параметр, тип целый);

JZ         - глобальный номер столбца матрицы Z, который указывает на начало
               подматрицы  sub (Z) (глобальный входной параметр, тип целый);

DESCZ - дескриптор распределенной матрицы Z (одномерный массив длины DLEN);
               если тип дескриптора одномерный ( DTYPE_Z = 501), то DLEN ≥ 7; если
               тип дескриптора двумерный ( DTYPE_Z = 1), то DLEN ≥ 9;
           DESCZ содержит информацию о размещении матрицы Z в памяти; полное
           описание DESCZ см. в п. 2.3.2; заметим, что DESCZ (CTXT_) должен быть равен
           DESCA (CTXT_); (глобальный и локальный входной параметр, тип целый);

WORK  - одномерный рабочий массив двойной точности длины LWORK; на выходе в
                элементе WORK (1) возвращается необходимая длина рабочего пространства
                (локальное рабочее пространство, локальный выходной параметр);

LWORK - задаваемая длина рабочего пространства WORK;
              LWORK ≥ 5 * N + MAX ( 5 * NN , NP0 * MQ0 + 2 * NB * NB) +
              ICEIL( NEIG, NPROW * NPCOL) * NN, где
              NP0 = NUMROC ( NN, NB, 0, 0, NPROW) - число строк локальной части
             матрицы sub (A);
             MQ0 = NUMROC ( MAX(NEIG, NB, 2), NB, 0, 0, NPCOL) - число столбцов
             локальной части матрицы sub (A);
             NEIG - число собственных векторов;
             NB = DESCA (NB_);
              NN = MAX( N, NB, 2 )
             Минимальная требуемая величина LWORK может быть вычислена самой
             подпрограммой, если к ней обратиться со значением LWORK = -1; при этом
             подпрограмма не производит никаких других вычислений, а требуемое значение
             LWORK возвращается в элементе массива WORK (1);
             (локальный входной параметр, тип целый);
     IFAIL    - обеспечивает дополнительную информацию, когда INFO ≠  0; если 
                     MOD (INFO/16,2) ≠ 0, то IFAIL содержит порядок наименьшего минора
                     матрицы  B, который не является положительно определенным (глобальный 
                     выходной параметр, тип целый);
INFO    - целая переменная, диагностирующая результат работы подпрограммы
               (глобальный выходной параметр)
               = 0 - успешное завершение работы;
               < 0 - если i-й фактический параметр подпрограммы является массивом и его 
                        j-й элемент имеет недопустимое значение, тогда INFO = - ( i * 100 + j );
                       если i - ый фактический параметр является скаляром и имеет
                       недопустимое значение, тогда INFO = - i;
>0 - если MOD(INFO/8,2).NE.0, то подпрограмма PDSTEBZ не смогла
        вычислить собственные значения;
        если MOD(INFO/16,2).NE.0, то матрица  B не является положительно
        определенной, а параметр IFAIL содержит порядок наименьшего минора,
        который не является положительно определенным.

Версии 

     PSSYGV2, PCHEGV2, PZHEGV2  -  вычисление всех собственных значений в
              обобщенной проблеме собственных значений для вещественных симметричных
              (эрмитовых) матриц A и B для случаев вещественных данных одинарной
              точности, комплексных данных одинарной точности и комплексных данных
               двойной точности соответственно

Вызываемые подпрограммы 

       Здесь указаны только базовые подпрограммы (т.е. подпрограммы второго уровня), которые вызываются из целевой подпрограммы (т.е. подпрограммы первого уровня).

    PDSYNGST, PZHENGST  -  приведение обобщенной проблемы собственных значений
       к стандартной форме

    PDSYEV2, PZHEEV2  -  вычисление всех собственных значений и собственных
        векторов симметричной (эрмитовой) матрицы в линейной проблеме собственных
        значений

    PDPOTRF, PZPOTRF  -   приведение симметричной (эрмитовой) матрицы к
        трехдиагональной форме методом отражений

    PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой

Замечания по использованию 

     1.  В подпрограммах PSSYGV2, PCHEGV2 и PZHEGV2 параметры A, B, Z и WORK
     имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметр W – 
     имеет тип REAL, REAL и DOUBLE PRECISION соответственно

     2. Список параметров подпрограммы PZHEGV2 имеет следующий вид:

PZHEGV2 (IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, W, Z, IZ, JZ, 
                    DESCZ, WORK, LWORK, RWORK, LRWORK, IFAIL, INFO)
 где дополнительные параметры RWORK и LRWORK означают следующее:

RWORK   - одномерный рабочий массив двойной точности длины LRWORK;
                  на выходе в элементе RWORK(1) возвращается необходимая длина рабочего
                  массива RWORK (локальное рабочее пространство, локальный выходной
                  параметр);

    LRWORK - задаваемая длина рабочего пространства RWORK;
                        LRWORK ≥ 4*N + MAX( 5*NN, NP0 * MQ0 ) +
                                                 + ICEIL( NEIG, NPROW*NPCOL)*NN , где
                   NEIG = число требующихся собственных векторов
                   NB = DESCA( NB_ )
                   NN = MAX( N, NB, 2 )
                   NP0 = NUMROC( NN, NB, 0, 0, NPROW )
                   MQ0 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL )
                   (локальный входной параметр, тип целый);

    3. Используются подпрограммы BLACS_GRIDINFO  из пакета BLACS,
     DSCAL из пакета BLAS,
     LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK1MAT, PCHK2MAT, PXERBLA из
     библиотеки ScaLAPACK_TOOLS

Пример использования

Вычисление всех собственных значений и собственных векторов в обобщенной проблеме собственных значений (1) для симметричных матриц A и B, которые имеют  следующий вид ( заданы только верхние треугольники ):

|   10    2    3    1    1    |
|    0   12    1    2    1    |
|    0    0   11    1   -1    |
|    0    0    0    9     1    |
|    0    0    0    0   15    |

|   12    1   -1    2    1    |
|    0   14    1   -1    1    |
|    0    0   16   -1    1    |
|    0    0    0   12   -1    |
|    0    0    0    0   11    |

     Порядок матриц N = 5. Матрицы разбиваются на блоки с размером блоков NB = 2. Осуществляется пропуск программы на четырех процессах, которые образуют решетку 2 на 2. Подробный разбор данного примера см. п. 3.4.

Фрагмент фортранного текста вызывающей программы
Полный текст этого теста приведен в Приложении 2, а принятые обозначения приведены  в п. 3.3.
Матрицы A и B инициализируются посредством обращения к подпрограмме PDELSET из библиотеки ScaLAPACK_TOOLS.

          PROGRAM TDSYGV21
          include ‘mpif.h’
INTEGER                         LWORK, MAXN, LDA, MAXPROCS, NOUT
PARAMETER                  ( LDA = 100, LWORK = 146, MAXN = 100,
               $                                         MAXPROCS = 512, NOUT =6 )
CHARACTER                 UPLO
PARAMETER                 ( UPLO = ‘U’ )
INTEGER                        CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB, IZ, JZ,
               $                                          NPCOL, NPROCS, NPROW, IA, JA, IB, JB, IBTYPE, IFAIL

INTEGER                        DESCA( 9 ), DESCB(9), DESCZ(9)
DOUBLE PRECISION   A( LDA, LDA ),  B( LDA, LDA ),  Z( LDA, LDA ), W( MAXN ),
               $                                        WORK( LWORK )

EXTERNAL                    BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
               $                                           BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
               $                                           BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
               $                                           PDSYGV2

        N = 5
        NB = 2
        NPROW = 2
        NPCOL = 2
        IA = 1
        JA = 1
        IB = 1
        JB = 1
        IZ = 1
        JZ = 1
        IBTYPE = 1

        CALL  BLACS_PINFO( IAM, NPROCS )
        IF( ( NPROCS .LT. 1 ) ) THEN
           CALL  BLACS_SETUP( IAM, NPROW*NPCOL )
        END IF

        CALL  BLACS_GET( -1, 0, CTXT )
        CALL  BLACS_GRIDINIT( CTXT, ‘R’, NPROW, NPCOL )
        CALL  BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )

         IF( MYROW .EQ. -1 ) GO TO 20

        CALL  DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
        CALL  DESCINIT( DESCB, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
        CALL  DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
* Построение матриц A и B

        CALL  PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
* Решение обобщенной проблемы собственных значений

        CALL  PDSYGV2( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, W, Z,
                                         IZ, JZ, DESCZ, WORK, LWORK, IFAIL, INFO )

        CALL  BLACS_GRIDEXIT( CTXT )
     20       CONTINUE
*
         CALL  BLACS_EXIT( 0 )
         STOP
         END

         SUBROUTINE  PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )

* PDMATINIT генерирует и распределяет матрицы A и B по решетке процессов
*
INTEGER                        IA, INFO, JA, N
INTEGER                        DESCA( * ), DESCB( * )
DOUBLE PRECISION   A( * ), B( * )
INTEGER                        I, J, MYCOL, MYROW, NPCOL, NPROW
EXTERNAL                    PDELSET

          INFO = 0

        IF( IA .NE. 1 ) THEN
               INFO = -3
        ELSE IF( JA .NE. 1 ) THEN
               INFO = -4
        END IF

         DO 2  J = 1, N 
         DO 1  I = 1, N 
               CALL  PDELSET( A, I, J, DESCA, 0.0D0)
               CALL  PDELSET( B, I, J, DESCB, 0.0D0)
1       CONTINUE
2       CONTINUE

         DO 3  I = 1, N-1 
                J = I + 1
                     CALL  PDELSET( A, I, J, DESCA, 1.0D0)
                   CALL  PDELSET( B, I, J, DESCB, 1.0D0)
            END IF
3       CONTINUE

          J = N
          DO 5  I = 1, N-1
                 CALL  PDELSET( A, I, J, DESCA, 1.0D0)
                 CALL  PDELSET( B, I, J, DESCB, 1.0D0)
5        CONTINUE

          CALL  PDELSET( A, 1, 1, DESCA, 10.0D0)
          CALL  PDELSET( A, 1, 2, DESCA, 2.0D0)
          CALL  PDELSET( A, 1, 3, DESCA, 3.0D0)
          CALL  PDELSET( A, 1, 4, DESCA, 1.0D0)
          CALL  PDELSET( A, 2, 2, DESCA, 12.0D0)
          CALL  PDELSET( A, 2, 4, DESCA, 2.0D0)
          CALL  PDELSET( A, 3, 3, DESCA, 11.0D0)
          CALL  PDELSET( A, 3, 5, DESCA, -1.0D0)
          CALL  PDELSET( A, 4, 4, DESCA, 9.0D0)
          CALL  PDELSET( A, 5, 5, DESCA, 15.0D0)

          CALL  PDELSET( B, 1, 1, DESCA, 12.0D0)
          CALL  PDELSET( B, 1, 3, DESCA, -1.0D0)
          CALL  PDELSET( B, 1, 4, DESCA, 2.0D0)
          CALL  PDELSET( B, 2, 2, DESCA, 14.0D0)
          CALL  PDELSET( B, 2, 4, DESCA, -1.0D0)
          CALL  PDELSET( B, 3, 3, DESCA, 16.0D0)
          CALL  PDELSET( B, 3, 4, DESCA, -1.0D0)
          CALL  PDELSET( B, 4, 4, DESCA, 12.0D0)
          CALL  PDELSET( B, 4, 5, DESCA, -1.0D0)
          CALL  PDELSET( B, 5, 5, DESCA, 11.0D0)

          RETURN
          END
Результаты:
Значение  INFO  =   0
Собственные значения:
W(1) =  0.432787211016963
W(2) =  0.663662748392315
W(3) =  0.943859004668386
W(4) =  1.10928454001752
W(5) =  1.49235323254300

Собственные векторы:
Z(1, 1) =  0.134590573961348964D+00
Z(2, 1) = -0.612947224715858080D-01
Z(3, 1) = -0.157902562211254954D+00
Z(4, 1) =  0.109465787723999147D+00
Z(5, 1) = -0.414730117966489037D-01

Z(1, 2) =  0.829198064861523532D-01
Z(2, 2) =  0.153148395665945264D+00
Z(3, 2) = -0.118603667911389210D+00
Z(4, 2) = -0.182813041785799935D+00
Z(5, 2) =  0.356172036818147315D-02

Z(1, 3) =  0.191710031573882250D+00
Z(2, 3) = -0.158991211513700409D+00
Z(3, 3) =  0.748390709386720365D-01
Z(4, 3) = -0.137468929466840417D+00
Z(5, 3) =  0.889778923495063573D-01

Z(1, 4) =  0.142011959884590616D+00
Z(2, 4) =  0.142419950546665203D+00
Z(3, 4) =  0.120997623004495347D+00
Z(4, 4) =  0.125531015188065975D+00
Z(5, 4) =  0.769220728307276440D-02

Z(1, 5) = -0.763867178777808409D-01
Z(2, 5) =  0.170980018713419862D-01
Z(3, 5) = -0.666645336709070502D-01
Z(4, 5) =  0.860480093056616990D-01
Z(5, 5) =  0.289433414168931813D+00

Приложение 2

Пример составления головной программы
для использования целевой подпрограммы из Приложения 1

          PROGRAM TDSYGV21
*        
*       Тест к подпрограмме   PDSYGV2
*
           include 'mpif.h'
*        Именованные константы - параметры задачи ..

INTEGER            LWORK, MAXN, LDA, MAXPROCS, NOUT
PARAMETER          (LDA = 100, LWORK=146, MAXN = 100,
*                     PARAMETER      ( LDA = 100, LWORK = -1, MAXN = 100,
                $                                    MAXPROCS = 512, NOUT =6)
CHARACTER          UPLO
PARAMETER          ( UPLO = ‘U’ )
*
*       Локальные переменные ..
*
INTEGER            CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB,
               $                              NPCOL, NPROCS, NPROW, IA, JA, IB, JB, IBTYPE
               $                              IZ, JZ, IFAIL
*
*        Локальные массивы ..
*
INTEGER                       DESCA( 9 ), DESCB( 9 ), DESCZ( 9 )
DOUBLE PRECISION   A( LDA, LDA ), B( LDA, LDA ), Z( LDA, LDA ),
              $                                          W( MAXN ), WORK( LWORK ), WORK1( 5 )
    $                                          W( MAXN ), WORK( 146 ), WORK1( 5 )
*
*        Внешние подпрограммы ..
*
EXTERNAL           BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
             $                                   BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
             $                                   BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
             $                                  PDSYGV2
*
*        Выполняемые операторы
*
*        Постановка задачи
*
      N = 5
     NB = 2
     NPROW = 2
     NPCOL = 2
     IA = 1
     JA = 1
     IB = 1
     JB = 1
     IZ = 1
     JZ = 1
     IBTYPE = 1
*
*       Инициализация пакета BLACS
*
     CALL BLACS_PINFO( IAM, NPROCS )
     IF( ( NPROCS.LT.1 ) ) THEN
          CALL BLACS_SETUP( IAM, NPROW*NPCOL )
     END IF
*
*       Инициализация контекста BLACS'а и решетки процессов
*
     CALL BLACS_GET( -1, 0, CTXT )
     CALL BLACS_GRIDINIT( CTXT, ‘R’, NPROW, NPCOL )
     CALL BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*        Если процесс не является частью данного контекста, то
*        переход на конец программы
*
     IF( MYROW.EQ.-1 ) GO TO 20
*
*        Инициализация дескрипторов для матриц  A, B и Z
*
     CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
     CALL DESCINIT( DESCB, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
     CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
*       Построение матриц  A и B
*
     CALL PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
*       Печать входных данных
*
     IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN 
          WRITE( NOUT, FMT = 92 )
         WRITE( NOUT, FMT = 98 )N, N, NB, NB
         WRITE( NOUT, FMT = 97 )NPROW*NPCOL, NPROW, NPCOL
         WRITE( NOUT, FMT = * ) ' Исходные данные '
         WRITE( NOUT, FMT = 88 ) UPLO
         WRITE( NOUT, FMT = * ) ‘  DESCA’
         WRITE( NOUT, FMT = 85 )DESCA
         WRITE( NOUT, FMT = * ) ‘  DESCB’
         WRITE( NOUT, FMT = 85 )DESCB
         WRITE( NOUT, FMT = * ) ‘ Матрица A’
     END IF
    85      FORMAT( / 9I5 /)
    88     FORMAT( / ‘ UPLO = ‘, A /)
*
     CALL PDLAPRNT( N, N, A, IA, JA, DESCA, 0, 0, ‘A’, NOUT, WORK1 )
*
     IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
          WRITE( NOUT, FMT = * ) ‘ Матрица B’
     END IF
     CALL PDLAPRNT( N, N, B, IB, JB, DESCB, 0, 0, ‘B’, NOUT, WORK1 )
*
*        Решение обобщенной проблемы собственных значений
*
     CALL PDSYGV2( IBTYPE, UPLO, N, A, IA, JA, DESCA,
     $                                       B, IB, JB, DESCB, W, Z, IZ, JZ, DESCZ,
     $                                     WORK, LWORK, IFAIL, INFO )
*
*       Печать результатов
*
     IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
          WRITE( NOUT, FMT = * ) ‘ Результаты ‘
          WRITE( NOUT, FMT = 96 )INFO
*                     WRITE( NOUT, FMT = 86 )WORK(1)
           WRITE( NOUT, FMT = * ) ‘ собственные числа’
           DO 10 I = 1, N
                  PRINT *, ‘ W(‘, I, ‘)=’, W( I ), ‘;’
  10        CONTINUE
          WRITE( NOUT, FMT = * ) ‘ собственные векторы’
     END IF
*
*        Печать собственных векторов
*
     CALL PDLAPRNT( N, N, Z, IZ, JZ, DESCZ, 0, 0, ‘Z’, NOUT, WORK1 )
*
     CALL BLACS_GRIDEXIT( CTXT )
  20     CONTINUE
*
     CALL BLACS_EXIT( 0 )
*
  99      FORMAT( ‘W=diag([’, 4D16.12, //’]);’ )
*
  92      FORMAT( / 'Тест к подпрограмме PDSYGV2' )
  98      FORMAT( / ' Вычисление всех собственных значений и',
     $                          ' собственных векторов '/
     $                          ' в обобщенной проблеме собственных значений ',
     $                          ' для матриц  A и B,'/
     $                          ' где A и B - вещественные симметричные матрицы ',
     $                            I2, ' на ', I2,',',
     $                          ' разбитые на блоки', I2, ' на ', I2, /)
*
  97        FORMAT( / ' Запуск на ', I2, ' процессах,',
        $                            ' образующих решетку размером ', I2, ' на ', I2 /)
  96        FORMAT( / ‘ Значение INFO  = ‘, I6 /)
* 86        FORMAT( / ‘ WORK(1) =’, D16.5 /)
     STOP
     END
*
     SUBROUTINE PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, 
       $                                                       INFO )
*
*           PDMATINIT генерирует и распределяет матрицы  A и B  по решетке процессов
*
*           Именованные константы – параметры
             INTEGER            CTXT_
             PARAMETER      ( CTXT_ = 2 )

*           Параметры - переменные
     INTEGER            IA, INFO, JA, N
*           Параметры - массивы
             INTEGER                        DESCA( * ), DESCB( * )
             DOUBLE PRECISION   A( * ), B( * )
*           Локальные переменные
     INTEGER            I, J, MYCOL, MYROW, NPCOL, NPROW
*           Внешние подпрограммы
     EXTERNAL           PDELSET
*           Внутренние функции
     INTRINSIC          DBLE
*
*          Выполняемые операторы
*
INFO = 0
*
IF( IA.NE.1 ) THEN
      INFO = -3
ELSE IF( JA.NE.1 ) THEN
      INFO = -4
END IF
*
DO 2 J = 1, N 
DO 1 I = 1, N 
     CALL PDELSET( A, I, J, DESCA, 0.0D0)
    CALL PDELSET( B, I, J, DESCB, 0.0D0)
1       CONTINUE
2      CONTINUE
DO 3 I = 1, N-1
       J = I+1
      CALL PDELSET( A, I, J, DESCA, 1.0D0)
      CALL PDELSET( B, I, J, DESCB, 1.0D0)
  3     CONTINUE
J = N
DO 5 I = 1, N-1
      CALL PDELSET( A, I, J, DESCA, 1.0D0 )
      CALL PDELSET( B, I, J, DESCB, 1.0D0)
  5      CONTINUE
*
CALL PDELSET( A, 1, 1, DESCA, 10.0D0)
CALL PDELSET( A, 1, 2, DESCA, 2.0D0)
CALL PDELSET( A, 1, 3, DESCA, 3.0D0)
CALL PDELSET( A, 1, 4, DESCA, 1.0D0)
CALL PDELSET( A, 2, 2, DESCA, 12.0D0)
CALL PDELSET( A, 2, 4, DESCA, 2.0D0)
CALL PDELSET( A, 3, 3, DESCA, 11.0D0)
CALL PDELSET( A, 3, 5, DESCA, -1.0D0)
CALL PDELSET( A, 4, 4, DESCA, 9.0D0)
CALL PDELSET( A, 5, 5, DESCA, 15.0D0)
*
CALL PDELSET( B, 1, 1, DESCB, 12.0D0)
CALL PDELSET( B, 1, 3, DESCB, -1.0D0)
CALL PDELSET( B, 1, 4, DESCB, 2.0D0)
CALL PDELSET( B, 2, 2, DESCB, 14.0D0)
CALL PDELSET( B, 2, 4, DESCB, -1.0D0)
CALL PDELSET( B, 3, 3, DESCB, 16.0D0)
CALL PDELSET( B, 3, 4, DESCB, -1.0D0)
CALL PDELSET( B, 4, 4, DESCB, 12.0D0)
CALL PDELSET( B, 4, 5, DESCB, -1.0D0)
CALL PDELSET( B, 5, 5, DESCB, 11.0D0)
*
RETURN
*
*       End of  PDMATINIT
*
END

ЛИТЕРАТУРА

1. Антонов А.С. Параллельное программирование с использованием технологии MPI:
Учебное пособие. М.: Изд-во МГУ, 2004.-71с.
2. Антонов А.С. Практический курс MPI. MPI: A Message-Passing Interface Standard (Version 1.1).
3. Библиотека численного анализа НИВЦ МГУ, http://num-anal.srcc.msu.ru/ 

4.Воеводин В.В. Математические модели и методы в параллельных процессах. -М.: Наука, 1986.-296с.

5. Воеводин В.В. Математические основы параллельных вычислений.  М.: Изд-во МГУ, 1991.-345с.
6. Воеводин В.В. Параллельные структуры алгоритмов и программ. М.: ОВМ АН СССР, 1987.-148с.
7. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления.  С.-П.-«БХВ-Петербург», 2002.-608с.

8. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления.  М.: Наука, 1984.-320с.
9. Воеводин В.В. Вычислительные основы линейной алгебры.  М.: Наука, 1977.-304с.  
10. Дацюк В.Н., Букатов А.А., Жегуло А.И. Введение в организацию и методы программирования многопроцессорных вычислительных систем (методическое пособие, часть I).  Ростов-на-Дону: Изд-во РГУ, 2000.-36с.
11. Дацюк В.Н., Букатов А.А., Жегуло А.И. Среда параллельного программирования MPI (методическое пособие, часть II).  Ростов-на-Дону: Изд-во РГУ, 2000.-65с.  
12. Немнюгин С.А., Стесик О.Л. «Параллельное программирование для многопроцессорных вычислительных систем». Петербург: БХВ-Петербург, 2002.-400с.  
13. Корнеев В.Д. Параллельное программирование в MPI. Новосибирск: Изд-во СО РАН, 2000.-213с.
14 E. Anderson, Z. Bai, C. Bischof, J.Demmel, J. Dongarra, J. Du Croz, A.  Greenbaum, S. Hammarling, A. Mckenney, S. Ostrouchov, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics: Philadelphia, 1995.
     15. L. S. Blackford, J. Choi, A. Cleary, J. Demmej, I. Dhillon, J. J. Dongarra, S. Hammarling, 
     G. Henry, A. Petitet, K. Stanley, D. W. Walker, and R. C. Whaley. ScaLAPACK: A portable
     linear algebra library for distributed memory computers: design issues and performance //
     Proceedings of Supercomputing ’96. New York: IEEE Press, 1996.
     http://www.supercomp.org/sc96/proceedings/).

     16. http://www.netlib.org/scalapack/index.html/.

17. J. Choi, J. Dongarra, and D. Walker. PB-BLAS: A Set of Parallel Block Basic Linear Algebra Subroutines // Practice and Experience, Vol. 8, 1996.  517-535.
18. http://www.netlib.org/scalapack/html/pblas_qref.html/.
19. J. Dongarra, R. van de Geun, and R.C. Whaley. Two dimensional basic linear algebra communication subprograms // Environments and Tools for Parallel Scientific Computing, Advances in Parallel Computing. Vol. 6. , Elsevier Science Publishers, 1993. Pp. 31-40. 
    20. http://www.netlib.org/blacs/ 
    21. G. Golub and C.F. Van Loan. Matrix Computations. Baltimore: Johns Hopkins University
    Press, 1996.
22. C.L. Lawson, R.J. Hanson, D. Kincaid, and F.T. Krogh. Basic linear algebra subprograms for Fortran usage // ACM Trans. Math. Soft.. Vol. 5. 1979. 308-323.
23. http://www.netlib.org/blas/ 
24. W. Lichtenstein and S.L. Johnsson. Block-cyclic dense linear algebra // SIAM J.  Sci. Stat. Comput. Vol. 14. 1993. 1259-1288.
25. L. Prylli and B. Tourancheau. Efficient block cyclic data redistribution // Proc. EUROPAR’96. Vol. 1. Lecture Notes in Computer Science. Berlin: Springer, 1996. 155-165.
26. Fernando K.V. and Parlett B.N. Accurate singular values and differential qd algorithms // Numer.Math. 1994. Vol. 67, No. 2. 191-230.
27. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980.